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Numerical analysis of influence of selected 
elements on effectiveness of streamline rudder 


Tomasz Abramowski, Ph. D. 
Tadeusz Szelangiewicz, Prof. 
West Pomeranian University of Technology, Szczecin 


ABSTRACT 


During designing steering gear for large fast transport ships (e.g. container carriers), shipowners usually 
put forward strong demands concerning ship manoeuvrability. It means that streamline rudders should 
be characterized by a high effectiveness, i.e. fast increasing values of lifting force in function of rudder 
angle and large values of lifting force related to rudder area. As gabarites of streamline rudder depend on 
а form and draught of stern part of ship s hull, an improvement of rudder effectiveness can be reached by 
an appropriate selection of rudder profile and application of additional elements to rudder blade. 
This paper presents results of numerical investigations (by using CFD methods) of hydrodynamic forces 
acting on rudder blades of the same gabarites but based on different profiles. Such calculations were also 
performed for selected rudder blades fitted with additional elements intended for the improving of rudder 
effectiveness. 


Keywords: streamline rudder; improvement of rudder effectiveness; computational fluid dynamics (CFD) 


INTRODUCTION 


During designing the streamline rudder its designer tends 
to make its effectiveness as large as possible. The effect can 
be reached by selecting an appropriate geometry of rudder, 
including its profile, as well as by applying additional elements 
which on the one hand are able e.g. to increase lifting force, 
and on the other hand e.g. to reduce tip losses. 

In the subject-matter literature can be found dimensionless 
results of model tests of hydrodynamic characteristics of 
airfoils or other profiles used in designing a given rudder of 
determined dimensions. However similar results of model tests 
of the rudders fitted with additional elements can be not always 
recalculated to a designed rudder of determined dimensions. 
Determination of effectiveness of a designed rudder can be 
hence performed: 

* either by conducting its model tests in a model basin, 
* or by using computational fluid dynamics methods 

(CFD). 

Model tests of rudders are much more expensive than 
numerical analyses which can be presently conducted both for 
models and full scale objects in view of rather easy access to 
high-performance computers, today. 


SCOPE OF THE NUMERICAL 
INVESTIGATIONS 


The numerical investigations were conducted in two 
phases: 
I. In the first phase hydrodynamic characteristics of the 
profiles (Fig. 1) applicable to designing the rudders, were 
calculated. Dimensions of the rudders and their blade area 


were the same as of the final rudder installed on B573 ship 
[1], which served as the initial object in the analyses in 
question. 


NACA HSVA HSVA IFS IFS IFS 
64,018 MP 73-20 MP 71-20 58-TR15 61-TR25 62-TR25 


Fig. 1. The profiles used for the investigated rudders 


II. In the second phase certain modifications were introduced 
to the selected rudders of the above given profiles, and 
then successive calculations of their hydrodynamic 
characteristics were made, namely: 

* for the rudder of NACA 0018 profile, fitted with 
a pivoting flap in the rear part of the profile, containing 
its trailing edge (Fig. 3a). In compliance with the subject- 
matter literature it was assumed that the deflection angle 
of the flap relative to rudder axis is equal to that of the 
main part of rudder blade (if the rudder angle was equal 
to 30° then the whole deflection angle of the flap was 
equal to 60°). 
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e for the rudder of IFS58-TR15 profile with an additional THE METHOD AND DOMAIN OF 
plate fixed at rudder top to reduce generation of tip COMPUTATIONS 
vortices (Schilling rudder, Fig. 3b). 


The computational model of the rudder prepared for 
numerical calculations is presented in Fig. 4. 


Fig. 4. The computational model of the rudder 


Fig. 2. The final rudder of B573 ship [1], without any additional elements 
of the area A, = 39.48 т? and the aspect ratio А = 1.79 . : . 
a) The computational domain (Fig. 5) was so arranged as to 
PS make it possible to assume different inlet angles of waterflow 


velocity to the rudder. The upper and lower surface limiting 
the domain was defined to be the plane of symmetry. The rear 
surface was defined to be the pressure outlet and the front 
surface (a half of cylinder) and the side ones constituted the 
inlet for waterflow velocity to the domain. 


780 


Fig. 5. The computational domain of rudders 


1185 


Fig. 3. Modification variants of the selected rudders: 
a) NACA 0018 rudder fitted with a pivoting flap, 
b) IFS 58-TR15 rudder fitted with an additional plate (Schilling s rudder) 


The modified rudders maintained the same contour and 


. . ribs Fig. 6. The computational domain mesh. Notation: the blue area — velocity 
dimensions as those of the B573 ship's rudder. inlet, the yellow one — plane of symmetry, the red one — pressure outlet 
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Fig. 7. The structural mesh in the boundary layer (view from the top of the 
rudder) Notation: the white area — the surfaces defining the rudder blade. 


For the domain a polyhedral mesh (Fig. 6 and 7) refined in 
the vicinity of the rudder and fitted with structural boundary 
layer on surfaces modelling the rudder, was used. 


CALCULATIONS FOR THE RUDDERS 
WITHOUT ANY ADDITIONAL ELEMENTS 


The calculations in question were conducted with the use 
of Fluent system. The example pressure distributions and 
streamlines for the rudder of IFS58-TR 15 profile are presented 
in Fig. 8, 9 and 10. 

In Fig. 11 are presented results of the calculations in the 
form ofthe lifting force L and drag force D versus rudder angle 


2.10e+03 
1.05e403 
0.00е+00 


for the rudders based on ће profiles shown in Fi 0.1 Fig. 8. The pressure distribution: a) and streamlines, b) for 0° angle of 


attack (the rudder of IFS58-TR15 profile) 
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Fig. 9. The pressure distribution and streamlines for 40° angle of attack (the rudder of IFS58-TR15 profile) 
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Fig. 10. The tip vortex generated on the rudder of IFS58-TR15 profile at 40° angle of attack. 


There are depicted streamlines, pressure profiles on iso-surface behind the rudder as well as velocity vectors on the surface 


CHOC)» CID a а а aaa A INIT 


Fig. 12. The pressure distribution and streamlines for the IFS58-TR15 
profile rudder at 40? angle of attack 
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Fig. 11. Values of the lifting force L and drag force D versus rudder angle 
for the considered rudders 


CALCULATIONS FOR THE RUDDERS 
FITTED WITH ADDITIONAL ELEMENTS 


Dimensions of the rudders with additional elements were 
the same as for the rudders without them. The domain for 
computation of hydrodynamic forces acting on the modified 
rudders as well as the scope of the calculations was the same 
as in the case of the rudders without any modification. The 
pressure distributions and streamlines for the Schilling rudder 
of IFS58-TR15 profile are presented in Fig. 12. 

The calculation results of the lifting force L and drag force 
D are presented in Fig. 13 for the NACA 0018 profile rudder 
fitted with the additional flap, and in Fig. 14 — for the Schilling 
rudder of IFS58-TR15 profile. 


12000 
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Fig. 13. The hydrodynamic characteristics of the rudder of NACA0018 
profile. Notation: І, D — for rudder with flap; L0, DO — for standard rudder 
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Fig. 14. The hydrodynamic characteristics of the rudder of IFS58-TR15 
profile. Notation: І, D — for Schilling rudder; L0,D0 — for standard rudder 


1. 


CONCLUSIONS 


The performed numerical experiment was aimed at 
demonstration of high usefulness of CFD methods in 
investigating and designing the streamline rudders. The 
calculations performed for rudders of different profiles 
showed how far selection of a profile is important for 
a designed ship to which various requirements concerning 
its manoeuvrability or limitations imposed on its rudder, 
may be assigned. The achieved results can be assessed only 
qualitatively (by comparing hydrodynamic characteristics 
for particular profiles) because results of model tests of 
the rudder installed on B 573 ship were lacking. There are 
available results of model tests of foils of given profiles 
but having a very large aspect ratio. However in the case 
in question the rudders of strictly determined dimensions 
have been analyzed. 

The two modified rudders analyzed in the second phase of 
the numerical experiment revealed much more favourable 
characteristics in contrast to the same rudders without any 
modification — however the obtained results should be 
considered rather qualitative, but not quantitive ones as 
additional investigations should be performed to obtain 
more exact results. In spite of that the quantitive changes 
have been found in line with expectations. 


3. 


The performed analysis showed that the NACA0018 
profile rudder fitted with the flap has much more favourable 
characteristics at small values of rudder angle. Probably, 
at smaller values of flow velocity in the whole domain and 
at its large values in the surounding of the rudder behind 
the propeller, the rudder has curved the flow behind the 
propeller changing direction of thrust, in consequence. Such 
situation can happen during ship’s manoeuvers at a very 
small (or even zero) speed of the ship. 

The Schilling rudder, i.e. the modified rudder of IFS 58- 
TR15 profile, exhibits also more favourable characteristics 
mainly at small values of angle of attack. The difference can 
be also observed at the angle values close to separation of 
flow where run of the characteristics is more stable than in 
the case of the standard rudder. This results from a reduced 
tip vortex effect; however, as observed in Fig. 12, the tip 
vortex has not been fully removed. 

In the drawings of this paper only certain examples of 
the performed computations of pressure distribution and 
streamlines, have been presented. The complete set of 
results of numerical analyzes in question can be found in 
the final report on the R&D project [2]. 
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A fatigue life calculation method for structural 
elements made of D16CzATW aluminium alloy 


Jozef Szala, Prof. 
Grzegorz Szala, Ph. D. 
University of Technology and Life Science, Bydgoszcz 


ABSTRACT 


In the calculating of fatigue life of structural elements the methods based on hypotheses of fatigue damage 
summation are commonly applied. In using the methods, apart from choice of an appropriate hypothesis, 
it is necessary to know a loading spectrum which usually constitutes a set of sinusoidal cycles of variable 
parameters (S , S), as well as fatigue characteristics in the form of Wöhler fatigue diagram, as a rule. 
Knowledge of the above mentioned data is low during structural design process as designed material 
objects and possibility of performing measurements are lacking. 
Hence in this work a calculation method of fatigue life of structural elements, based on the relations between 
fatigue life diagrams determined in constant-amplitude conditions (Wöhler diagrams) and those obtained 
under programmed or random loads, is described. The proposed method was experimentally verified by 
using results of fatigue tests on structural elements made of DI6CzATW aluminium alloy. 


Keywords: fatigue life, structutral elements, DIGCzATW aluminium alloy 


INTRODUCTION 


Lack of a coherent theory of fatigue process of materials 
including metals and their alloys makes that for engineering 
applications are used phenomenological hypotheses on 
summation of damages, whose description can be found 
a.o. in the monographs [1] and [2]. Out of many hypotheses, 
the oldest hypothesis, most comprehensively verified and 
commonly used is that based on linear summation of fatigue 
damages, given by Palmgren and Miner (P-M). In many cases 
its experimental verification showed a low conformity between 
results of fatigue life calculations and experimental tests, that 
has resulted in undertaking further efforts to elaborate new 
calcultation methods [3, 4, 5, 6]. For the calculations with the 
use of the damage summation hypotheses it is necessary to 
know additionally operational loads (stresses or strains) given 
in the form of loading spectra, as well as an appropriate fatigue 
diagram. Such loading spectra are elaborated by applying the 
cycle - counting methods described in [7], in order to obtain 
a set of sinusoidal cycles equivalent, as regards fatigue life, to 
random loading. The third item necessary for the calculations 
is a Wóhler fatigue diagram of parameters appropriate to 
a type of loading spectrum. The diagram can be determined 
experimentally or, if not possible, selected from catalogues 
or literature sources [8, 9]. Out of the above mentioned main 
elements of the calculating algorithm, the largest uncertainty 
is associated with the fatigue damage summation hypothesis. 
Experimental verification of P-M hypothesis, performed by 
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means of a few hundred tests revealed the difference between 
results ofthe experiments and calculations ranging from 0,3 to 
30-fold [10, 11]. The drawbacks of the method have been very 
clearly felt in structural design process where accurate data 
on operational loads, material cyclic properties and damage 
summation hypothesis, are lacking. 

It is then necessary to elaborate, excepting damage 
summation hypothesis, a simple calculation method based on 
assessment of main loading spectrum parameters and basic 
cyclic properties of material. To the above mentioned conditions 
corresponds the method based on statistical relations between 
fatigue life determined under sinusoidal load conditions 
(described by Wóhler diagram) and that determined under 
random load conditions or programmed ones (described by 
Gassner diagram). 

As observed in the fatigue life tests of structural elements 
under programmed load conditions (block spectra, usually) 
corresponding to operational loading, their test results depended 
on a level of variable loads and shape of loading spectra. Such 
relationship for steel elements was demonstrated in [10] and 
[12], and for aluminium elements - in [13]. The problem was 
also described in the monograph of Haibach [14] and the 
information bulletins of Germanischer Lloyd Group, a ship 
classification society, [15] and [16]. 

As proved in the experimental verification ofthe calculation 
method based on statistical relations between Wóhler and 
Gassner diagrams, performed on steel specimens and structural 
elements, the calculation results obtained with its use, have been 


closer to results from experimental tests than those achieved 
from the calculations to which Palmgren-Miner hypothesis on 
damage summation was applied, [17]. 

In the work [17] was analyzed an influence of a shape 
of loading spectrum on results of experimental verification 
of calculation method of fatigue life of steel specimens and 
elements. 

This work is aimed at extending the method onto aluminium- 
alloy elements, in which data determined from plain specimens 
are applied to calculation of geometrically notched elements 
and structural elements. 


FORMULATION OF THE PROBLEM 


The method proposed in this paper is based on the following 
assumptions: 
a) damage cumulation process is qualitatively independent on 
loading spectrum, 
b) fatigue damage cumulation depends quantitatively on 
parameters of loading spectrum (its shape and maximum 
load), 
constant-amplitude loading is a special case of loading 
spectrum elaborated on the basis of random loads, e.g. block 
spectrum, 
d) if number of cycles in a loading spectrum (for a specimen 
or structural element) is reduced to one, then fatigue life 
diagrams intersect each other in the same point; it means 
that the spectral form loses its sense, 
for different loading spectra of the same maximum load 
value, fatigue life of a specimen or structural element 
depends mainly on a form of a spectrum. 


c 


— 


е 


— 


In Fig. 1 the sinusoidal loading, random loading and loading 
spectra corresponding to them are schematically presented. 

Loading spectrum form is characterized by the diagram 
filling factor calculated from the formula as follows [7]: 


1 k 
С ^g È Sat: (1) 


which, in relation to the schematic spectrum diagram, is 
equivalent to the ratio of the areas F, (under the block spectrum 
envelope) to the area of the rectangle having the sides: S... 
and n . As results from the scheme (b) in Fig. 1 and the formula 
(1), the spectrum filling factor G is equal to 1.0 for sinusoidal 
loading, and for random loading amounts to $ < 1.0. The 
smaller value of the diagram filling factor the slighter loading 
conditions for structural element; an example of mutual location 
of the fatigue-life diagrams for various values of the diagram 
filling factor is presented in Fig. 2. 

The schematic procedure of the calculation complying 
with the above specified assumptions is presented in Fig. 3. 
The Wöhler diagram (W) and Gassner diagram (С) which 
intersect each other in the point с = c,, are described by the 
expressions: 


amax i-l 


- for the Wöhler diagram: 
1 
logS, = -—logN+c, (2) 
m, 
- for ће Gassner diagram: 
| 
logS, max = =—log Ne +С (3) 
m 


Smax 


Fig. 1. Schematic loading diagrams and relevant loading spectra: a) sinusoidal loading, 
b) sinusoidal loading spectrum, c) random loading, d) random loading spectrum 
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Fig. 2. Schematic diagram of fatigue life dependence 
on loading spectrum form 


Nea No 


N, Ne 
Fig. 3. Schematic diagram for fatigue life calculation 


From the scheme of Fig. 3 results the possibility of 
determining Gassner diagram on the basis of knownledge of 
Wöhler fatigue diagram and an arbitrary point of the Gassner 
diagram under determination.The fact is of a great importance 
for shortening time of tests because it requires to perform 
programmed fatigue tests (which are time-consuming) for only 
one point of Gassner diagram (a few fatigue tests at only one, 
of maximum load, level of the spectrum). 

For Sima = S, and c, = с 1, by comparing the formulae (2) 
and (3) to each other the following is achieved: 


m 
log Ne = — 2 (4) 
0 


The relation m/m, in function of the diagram filling factor 
б, for different kinds of steel specimens and elements, was 
experimentally determined in the following form, [10, 17]: 
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m rd (5) 


The power exponent r is contained within the range of 0.2 
7 0.4. It should be stressed that if one ofthe lower values ofthe 
range is assumed then conservative results (1.e. lower fatigue 
life values) are obtained. 

On determination of N from the formula describing Wöhler 


diagram is obtained: M 
C 0 
к= е (6) 
S, 


and on substitution of the formula (5) to the formula (4) the 
following is obtained, after transformations: 


Co (7) 


amax 

It should be explained that c, which appears in the formula (2), 
is equal to log C, which appears in the formulae (6) and (7). 

The experimental verification of the desribed calculation 
method, performed for specimens and bonded structural joints 
for three values of the diagram filling factor © (0.34; 0.56 and 
0.77), showed a conformity of the results obtained from the 
calculation in question with those from the tests higher than 
in the case of the results from the calculation where selected 
fatigue damage summation hypotheses were applied. 


TESTS AND CALCULATION RESULTS 
FOR SPECIMENS AND STRUCTURAL 
ELEMENTS MADE OF D16CZATW 
ALUMINIUM ALLOY 


log No =с ‘m, log 


Method and program of the tests 


From the target of this work the necessity results to perform 
basic tests aimed at determination of static and cyclic properties 
of material used for objects to be tested, as well as verification 
tests of the calculation method of fatigue life of elements under 
variable random loads. 

In the first group of the tests, apart from the monotonic 
tension test, Wóhler fatigue diagrams have to be determined, 
as well as the test aimed at determination of value of the 
power exponent r appearing in the formula (5), has to be 
performed with the use of plain specimens under random 
loading conditions. 

The determination method of value of the power exponent 
r results from Fig. 3. For tests under random or programmed 
loading it is enough to have aWóhler diagram, to determine 
a point on Gassner fatigue-life diagram, the point A in Fig. 3, 
and, to form a complete Gassner's diagram of fatigue life by 
drawing the straight line containing the points A and c. 


Time 


Fig. 4. Random loading assumed for the tests, 
according to the publication [18] 


Variable loading class intervals 


Diagram of cumulated numbers 
of local minimum points 


0.1 1.0 


Zero-level of loads or stresses 


10 100 % 


Number of local extremes 


Fig. 5. Loading spectrum for the random loading of Fig. 4, obtained by using the method of counting local extremes 


Variable loading and its spectrum 


As an example operational loading the data contained in 
the publication FALSTAFF-130 [18] were taken. Its graphical 
form is shown in Fig. 4, and Fig. 5 contains its spectrum 
achieved by applying the local extremes counting method [7]. 
On the ordinate axis of the loading diagram, in compliance 
with the assumption of the FALSTAFF project, are given class 
intervals within which local extremes have been counted. Such 
approach to loading spectrum makes performing the tests on 
various levels of the maximum stress S „p easier, necessary 
to determine fatigue life diagram for a tested object. The tests 
were conducted on a few levels of the maximum stress 5 of 
the spectrum. The spectrum filling factor calculated from the 
formula (1) was equal to б = 0.375. 


The testing objects 


The testing objects made of D16CzATW aluminium alloy 
are shown in Fig. 6. Their basic mechanical properties including 
cyclic ones were determined with the use of the flat unnotched 
(plain) specimens whose shape and dimensions are given in 
Fig. 6a. The tests for verification of the calculation method in 
question were conducted on the plate specimens with rivet holes, 
Fig. 6b, as well as on the riveted structural element, Fig. 6c. The 


chemical composition of D16CzATW alloy of which the tested 
objects were manufactured, is given in Tab. 1. 


Results of the tests 


The monotonic tension test was conducted on the standard 
specimens acc. PN-EN 10002-1+AC1 standard (Fig. ба) with 
the use of Instron 8502 test machine equipped with suitable 
instruments and programs. The test results are presented in 
Tab. 2. 

In Tab. 3 are given the data concerning fatigue properties of 
D16CzATW aluminium alloy in the range of high-cyclic fatigue 
(HCF) under sinusoidal pulsating load conditions. 

The complete Wóhler's diagram for the plain specimens, 
together with the depicted points of experimental results, is 
presented in Fig. 7. 

The pulsating load conforms the best to the operational load 
assumed for the tests in question (Fig. 4). 


Experimental determination of value of the 
power exponent r 


In order to obtain a value of the power exponent r the 
experiment was performed in compliance with the method 
above described in p. 3.1. To this end, the tests were conducted 


Tab. 1. Chemical composition of DIGCzATW aluminium alloy 


Pure alloy (Cz). cladded with AD1 plating of normal thickness (A). in supersaturated and natural ageing state 
(T). intended for covering (W) 


Tab. 2. Static properties of plates made of DIGCzATW aluminium alloy 


Rm Re Rm E A Z 
[MPa] [MPa] [MPa] [MPa] [96] [v6] 


68402 


Tab. 3. Cyclic properties of plates made of DIGCzATW aluminium alloy 


Wohler equation Fatigue limit [MPa] | Limit number of cycles 


Pulsating ШИ m- 53107 _ E 
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a) 120 


Fig. 6. The tested objects: a) the plain specimen prepared acc. PN-74/H-04327 standard, b) the specimen with rivet holes, c) the riveted structural element 


2164.7 


log N = 5.3107 log 


~ 


max 


103 104 105 106 
Number of cycles N [-] 


Fig. 7. Wöhler diagram for DIGCzATW aluminium alloy, determined under 
pulsating stress conditions (R — 0) 
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on six specimens under random loading conditions (Fig. 4) 
on the mean level of the maximum stress S axa = 350 [MPa]. 
The fatigue-life test results are presented in Tab. 4. The mean 
fatigue life value amounts to N, — 0,953. 

By transforming the formulae (4) and (5) the following is 
obtained: 


За logN (8) 
log N 
hence: BT 
r= log logN _1_ (9) 
Іов № / logg 


For $, „a= 350 [MPa], m= 5.3107, and C= 2164.7 [MPa] 
(Tab.3), from the formula (6) N = 1.56:10* is calculated. 


Tab. 4. Results of the fatigue-life tests under random loading conditions, for the stress value S. 


No. of specimen 


= 350 [MPa] 


max 


Mean fatigue life value N, 


Fatigue life N., 


0.953 


As results from the data contained in p. 3, the spectrum filling 
factor amounts to G = 0.375. 

On substitution of the above mentioned data to the formula 
(9) г = 0.36 is obtained. 

The value is contained within the variability limits r=0.2 + 0.4 
determined for steel specimens and elements, [10]. 


Experimental verification of results 
of fatigue life calculation 


Fatigue life calculations and tests of plain specimens 


Comparison of results of calculations with those of fatigue 
life tests on the plain specimens made of DIGCzATW alloy 
is mainly aimed at verification of the calculation method 
assumptions formulated in p. 2 of this work, throughout the 
whole variability range of the stresses 5 in the loading 
spectrum. 

The following data were asssumed for the calculations: 

- the spectrum filling factor — in accordance with p. 3.2: 

670.375 
- the power exponent — in accordance with p. 3.5: r= 0.36 
- the Wohler diagram — in accordance with the data of Tab. 3 

- described by the formula: 


log N = 5.3107 log 2а 


(ба) 
тах 


The formula (7) for fatigue life calculation takes the form: 
log М?” 21.4235-log N (4a) 


Results of the calculations are presented in Tab. 5, row 4 
and 5. 

The tests of the plain specimens made of DIGCzATW alloy 
were peformed under the random load conditions described 
in p. 3.2. The tests were conducted on 30 specimens, for 9 
loading levels, namely: 450, 440, 400, 360, 350, 320, 300, 
280 and 250 [MPa]. 

Results of the tests together with the Gassner diagram based 
on them are presented in Fig. 8. 

The diagram is described by the formula: 


log МЕ = 7.4184 log — 


(10) 
тах 
Results of the calculations of N.* values from the formula 
(10) for selected levels of 5, are presented in Tab. 5, row 
6 and 7. 


200 ——— áÓ ———-—— , —————————— 


105 106 10 
Number of cycles N, [-] 


Fig. 8. Gassner fatigue life diagram (a) for plain specimens made 
of DIGCzATW aluminium alloy 


Comparison of the Gassner diagrams determined on the 
basis of the calculations and tests is shown in Fig. 9. 

From the diagrams of Fig. 9 result only slight differences 
between the calculated fatigue life and the experimentally 
determined. The relative deviation calculated from the formula 
(11) amounts — in extreme cases — to: à, = 5.4% for S „7450 
[MPa], б, = 10% for S „= 150MPa. 


max 
5 N 


ех obl 
_ Ne - № 


ех 


100% (11) 

The à, values which do Hot exceed 10% in relation to the 
results ofthe fatigue life calculations, should be considered very 
small. The very high calculation accuracy results mainly from 
the way of determining the power exponent r in the formula 
(8) on the basis of the experiment for the mean values of S... 
However this analysis is aimed at showing that in the whole 
5 ax Variability range from 150 to 450 [MPa] the assumptions 
for the calculation method are correct. 


Fatigue life calculations and tests of the notched 
specimens (acc. Fig. 6b) 


The analysis described in this point is mainly aimed 
at demonstrating that the data concerning the exponent r 
determined from standard plain specimens can be used as the 
basis for the calculations in the case of notched specimens and 
complex structural elements. 

The testing and calculating methods are the same as in p. 
3.6a, it means that Wóhler fatigue diagram, loading spectrum 


Tab. 5. Results of the fatigue life calculations and tests of the plain specimens made of DI6CzATW alloy 


150 


200 


300 350 


6.1443 


5.4823 


4.5490 4.1944 


1.4:106 


3.03610? 


3.54-10* 1.56:10* 


8.7427 


7.8041 


6.4755 5.9707 


5.53:10* 


6.37107 


2.99:]06 9.35`10° 


8.7001 


7.7759 


6.4728 5.9777 


5.014:10* 


59710 


2.87°10° 9.5°10° 
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Plain specimen, material: D16CZATW aluminium alloy 
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Fig. 9. Fatigue life diagrams for plain specimens made of D16CzATW aluminium alloy: a) experimentally determined under random loading conditions, 
b) based on calculations acc. formula (4a), shown against the Wohler diagram (c) 


parameters as well as values ofthe power exponent r constitute 
the basis for the calculations. The experimental tests of the 
elements were conducted under the random loading conditions 
described in p. 3.2. 

The Wóhler fatigue diagram for the plate with rivet holes, 
made of DIGCzATW aluminum alloy, is shown in Fig. 10. 


log N = 3.8095 log Ll 


max 


——— 


Piot 105 2-105 
Number of cycles N [-] 
Fig. 10. The Wöhler fatigue diagram for the plate with rivet holes 
(acc. Fig.6b), made of DIGCzATW aluminum alloy 


The diagram is described by the formula: 
3736 
log N 23.8095.1og 


(12) 


Results ofthe fatigue life tests under operational loading are 
shown in the Gassner diagram form in Fig. 11, and the formula 
(13) provides its description. 

Fig. 11. The Gassner fatigue diagram for the plate with rivet 
holes (acc. Fig.6b), made of DIGCzATW aluminum alloy 


log № = 4.0521og esc (13) 


max 

The results of the calculations and tests are collected in 
Tab. 6, like in the case of the plain specimens. The data given 
in the row 4 and 5, Tab. 6, were calculated on the basis of the 
formula (4a). 

The results of the calculations and tests given in Tab. 6 are 
graphically presented in the diagrams of Fig. 12. 

The maximum relative differences between the calculation 
and test results, determined from the formula (11) amount 
to: à, = -124% for S a = 150 [MPa], and à, = 42% for 


N 


$ = 400 [MPa]. 


тах 


Tab. 6. Results of the fatigue life calculations and tests of the plates with rivet holes (acc. Fig.6b), made of DIGCzATW alloy 


150 


250 300 


5.3193 


4.4741 4.1725 


2.9°10° 


2.98:10* 1.49:10* 


7.5720 


6.3689 5.9396 


3.7310’ 


2.34:10° 8.7:105 


7.2206 


6.3217 6.0008 


1.66:107 
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2.1:106 1.07106 


Specimen with holes, material: D16CZATW aluminium alloy 


Stress б, [ MPa] 


N 
Ре МЧА 
OCI Sot • РП 
ШП!! EET ТШ Г ШИ! 
БИЕ 


100 
1.E+03 1.E+04 1,E+05 1.E+06 1.Е+07 1,Е+08 1,Е+09 
Number of cycles М 


Fig. 12. Results of the fatigue life calculations and tests of the plates with rivet holes (acc. Fig. 6b), made of DIGCzATW aluminum alloy: 
а) experimentally determined under random load conditions, b) based on the calculations, shown against the Wöhler diagram (c) 


The above presented differences in the results of fatigue 
life analysis are relatively small as compared with the results 600 


available by means of other methods. 500 21. 21548 
log N= УВ 


Fatigue life calculations and tests 
of the riveted elements (acc. Fig. 6c) 


As the run of calculations and tests of the plates with rivet 
holes and the riveted elements is identical, in this point only 
the relevant formulae, diagrams and table are included. 

In Fig. 13 the Wóhler diagram for the riveted element acc. 
Fig. бс is presented, and in Fig. 14 — the results of the tests 
under operational loading conditions. 


The diagrams of Fig. 13 and 14 are respectively described 10° 104 10° MF 107 
as follows: Number of cycles N [-] 
- the Wóhler diagram Fig. 13. The Wöhler diagram for the riveted element 
2154.8 (acc. Fig. 6c) made of DIGCzATW aluminum alloy 
log N = 5.048 log (14) 
8 The fatigue life according to the proposed method was 
- the Gassner diagram calculated by using the formula (4a). 
6912 In the diagrams of Fig. 15 the results of the calculations 
log Ne = 4.529 log —— (15) апа tests given in Tab. 7 are graphically compared with the 
ee Wohler diagram. 


Tab. 7. Results of the fatigue life calculations and tests of the riveted elements (acc. Fig.6c), made of DIGCzATW alloy 
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400 


А 6912 
log М“ = 4509109 £ 


300 


Stress Sna | MPa] 


Е 5 6 7 7 „7 
5-10 10 10 2-0 3-10 
Number of cycles N, [-] 


Fig. 14. The Gassner fatigue life diagram for the riveted element 
(acc. Fig. бс) made of DIGCzATW aluminum alloy 


As results from the comparison, the greatest differences 
between the experimentally determined fatigue life and that 
calculated appear at the ends of the stress range used for the 
analysis in question. The relative values of the differences 
calculated according to the formula (11) amount to: à, = 0.37 
for S... — 250 [MPa], and à, = -5.08 for S „= 150 [MPa]. 

The calculated differences for the riveted element are greater 
than those for the plate with rivet holes. The fact results from 
the greater spread of results of fatigue life tests of the riveted 
element which is deemed incommensurably more complex 
from the point of view of run of fatigue process depending on 
various factors. 

If the spread band of the test results is taken into account 
on the confidence level of 0.95 then the fatigue life calculation 
results are contained within the band. 


SUMMARY 


The discussed calculation method for fatigue life of structural 
elements, based on statistical relations between the fatigue 
diagrams determined under constant - amplitude loading 
conditions (Wóhler diagrams) and those determined under 
random loading conditions characteristic for operational 
conditions of structural elements (Gassner diagrams), 
constitutes an essential contribution to the calculation 
methods based on phenomenological hypotheses of fatigue 
damage summation. 

Formulation of the assumptions for the proposed method 
is described in the publications [2] and [17] which also 
contain the experimental verification of the method, based 
on the comparison between the calculation results and those 
obtained from the tests on steel elements, performed in 
the loading conditions characterized by different values of 
the spectrum filling factor G (0.34; 0.56; 0.77). The crucial 
conclusion resulting from the analysis is that the analyzed 
calculation method may be effectively applied to load 
spectra irrespective of values of their factor С. 

In this work the applicability analysis of the specimens and 
structural elements made of the DI6CzATW aluminum alloy 
used in aircraft industry, has been performed. And, its original 
element, as compared with earlier publications, is the proof 
that the parameters which appear in the formula for calculating 
the fatigue life determined from standard specimens, are also 
applicable to calculations of complex structural elements; 
the statement mainly concerns value of the power exponent 
r appearing in the formulae (5) and (7). 

The comparative analysis ofthe calculation and experimental 
test results showed a satisfactory conformity and possible 
practical application of the proposed method to structural 


Riveted element, material: DIGCzATW aluminium alloy 
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Fig. 15. Graphical comparison of the results of the calculations and tests presented in Tab.7: a) experimentally determined under random load conditions; 
b) based on the calculations; shown against the Wöhler diagram (c) 
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design process when, as real material objects are lacking, 
only scarce data on operational conditions (data on 
measured loads) and on design features of such objects 
(being in statu nascendi), are at one’s disposal. 


This work has been elaborated in the frame of the research 
project No. 7 T07B 01018 financially supported by the Ministry 
of Science and Higher Education. 


NOMENCLATURE 


i — spectrum filling factor 
5 — stress (general notation), [MPa] 
S — maximum. stress: in a sinusoidal cycle, an 
analyzed section of random stresses (e.g. 
a measuring section of operational loading), 
loading spectrum, or loading programme, 
[MPa] 
— minimum stress: in a sinusoidal cycle, an 
analyzed section of random stresses (e.g. 
a measuring section of operational loading), 
loading spectrum, or loading programme, 
[MPa] 
S — mean stress value: in a sinusoidal cycle, an 
analyzed section of random stresses (e.g. 
a measuring section of operational loading), 
loading spectrum, or loading programme, 
[MPa] 
— stress amplitude in a sinusoidal cycle 
— stress amplitude in loading cycles: sinusoidal 
of i-th step, a loading spectrum or loading 
programme 
i — number ofa step in a loading spectrum or 
n n;i loading programme 
— fraction of number of S - amplitude cycles in 
total number of cycles appearing in a loading 
spectrum or programme 
n. = number of S, - amplitude cycles appearing in 
a loading spectrum 
n — total number of cycles within a loading 
Spectrum 
n.. = number of S, - amplitude cycles within 
a loading programme period 
n — total number of cycles within a loading 
programme 
À — number of repetitions of a loading programme 
period up to fatigue fracture 
— number of cycles (general notation) 
— number of cycles to fatigue fracture in 
sinusoidal S, — amplitude loading conditions 
N, — fatigue life 
F(S) — distribution of values of sinusoidal cycles 
m — parameter in the formulae which describe 
Wohler diagram 
Cim — parameters in the formulae which describe 
Gassner fatigue life diagram, Eq. (3) 
C — constant value in the formula which describe 
Wohler diagram of exponential form 
S — fatigue limit at the stress ratio R — 0 (pulsating 
load) 
ô — relative difference between fatigue life test and 
calculation results (see Eq. (11)). 


max 


min 


Indices 
ex — stands for operational results or values 
obl — stands for calculation results or values. 


BIBLIOGRAPHY 


1. Manson S.S.,Halford G.R.: Re-Examination of Cumulative 
Fatigue Damage Analysis — an Engineering Perspective. 
Engineering Fracture Mechanics, Vol.25, No. 5/6, 1986 

2. Szala J., Boroński D.: Material fatigue state assessment in 
diagnostics of machines and devices (in Polish). Wydawnictwo 
ITE-PIB (Publishing House of Institute of the Technology of the 
Exploitation - State Research Institute), Bydgoszcz, 2008 

3. Sevensson T., Johannesson P., de Mare J.: Fatigue life prediction 
based on variable amplitude tests-specific applications. 
International Journal of Fatigue, 27, 2005 

4. Liu Y., Mahadevans S.: Stochastic fatigue damage modeling 
under variable amplitude loading. International Journal of 
Fatigue, 29, 2007 

5. Buch A.: Prediction of the comparative fatigue performance 
for realistic loading distributions. Prog. Aerospace Sci. Vol. 33, 
1997 

6. Sonsino C.M.: Fatigue testing under variable amplitude 
loading. International Journal Fatigue, 29, 2007 

7. Kocanda S., Szala J.: Essentials of fatigue calculations (in 
Polish). PWN (State Scientific Publishing House), Warsaw, 1997 

8. Soprotivlenie ustalosti metalom i splavov. Kiev Naukowe 
Dumka, 1987 

9. Boyer H. E.: Atlas of Fatigue Curves. American Society for 
Metals (ASM), Ohio, 2003 

10.Szala J.: Fatigue life assessment of machines under random 
and programmed loading conditions (in Polish). Monograph, 
Mechanika 22, Wydawnictwo ATR (Publishing House of 
Technical Agricultural Academy), Bydgoszcz, 1980 

11.Lehman R.: Über die Gültigkeit der Hypothese der Linearen 
Schadens — akumulation. Maschinenbautechnik 5, 1970 

12.Haibach E.: Probleme der Betriebsfestigkeit von metalischen 
Konstruktionsteilen. VDI- Z133, 1971, Nr 5, 1.297/403 

13.Osterman H.: Die Lebensdauerabschdtzung bei 
Sonderkollektiven. LBF-Bericht Nr TB-80, 1968 

14.Haibach E.: Betriebs-Festigkeit. VDI Verlag GmbH, Düsseldorf, 
1989 

15.Fricke W., Petershagen, Paetzold H.: Fatigue Strength of 
Ship Structures, Part I: Basic Principles. GL-Technology, 
Information from Germanischer Lloyd Group, Nr 1/97, 
Hamburg, 1997 

16.Fricke W., Petershagen, Paetzold H.: Fatigue Strength of Ship 
Structures, Part II: Examples. GL-Technology, Information 
from Germanischer Lloyd Group, Nr 1/98, Hamburg,1998 

17.Szala J.: On a fatigue calculation method for structural elements 
under stochastic loading conditions (in Polish). Archiwum 
Budowy Maszyn (Archive of Machinery Engineering), Vol. 
XXIX, No. 3—4, 1982 

18.Szala J.: Selected fatigue-life problems of elements of ship 
machines and devices (in Polish). Proc. the Jubilee Scientific 
Conference: Research & Development — a chance for Polish 
shipbuilding industry, Ship Design and Research Centre, 
Gdansk, 2001. 


CONTACT WITH THE AUTHORS 


Józef Szala, Prof., 
Grzegorz Szala, Ph. D. 

Faculty of Mechanical Engineerig, 
University of Technology and Life Science, 
Prof. S. Kaliskiego 7 
85-763 Bydgoszcz, POLAND 
e-mail: jszpkm@utp.edu.pl 


POLISH MARITIME RESEARCH, No 3/2010 17 


POLISH MARITIME RESEARCH 3(66) 2010 Vol 17; pp. 18-25 
10.2478/v10012-010-0024-1 


Theoretical and mathematical models of the 
torque of mechanical losses in a hydraulic 
rotational motor for hydrostatic drive 


Zygmunt Paszota, Prof. 
Gdansk University of Technology 


Abstract 


The paper presents theoretical and mathematical models of the torque of mechanical losses 

in a hydraulic rotational motor with constant capacity qMt per one shaft revolution (with 

constant theoretical working volume V Mt) and with variable capacity qMgv = bM qMt per 

one shaft revolution (with variable geometrical working volume VMgv). The models are 

to be used in the laboratory and simulation investigations of motor energy losses aimed 
at evaluation of the motor energy efficiency and hydrostatic drive efficiency. 


Keywords: hydrostatic drive; hydraulic motor; energy efficiency 


INTRODUCTION 


The paper is a continuation of references [1+10], aimed 
at developing a method of evaluation of the losses and 
energy efficiency of the hydrostatic drives and displacement 
machines used in them. The method is based on theoretical and 
mathematical descriptions of losses in the pumps, hydraulic 
motors and in other elements of a hydrostatic drive system. 

Description of the hydraulic motor losses and energy 
efficiency is based on the diagram of the direction of increasing 
power stream in a rotational hydraulic motor, which is 
introduced instead of the Sankey diagram (Fig.1 [10]). 

The aim of the paper is to present theoretical and 
mathematical models of torque of mechanical losses in the 
rotational hydraulic motor ,,shaft — working chambers” 
assembly. The motor is a displacement machine with constant 
theoretical capacity q,,, per one shaft revolution (with constant 
motor theoretical working volume V,,,) or with variable 
geometrical capacity duas buau per one shaft revolution (with 
variable motor geometrical working volume V,, ww)" 

The models are to be used in the laboratory and simulation 
investigation of the motor energy losses, motor energy 
efficiency and hydrostatic drive efficiency. 


THEORETICAL MODELS OF THE 
TORQUE Му OF MECHANICAL LOSSES 
IN THE MOTOR ,SHAFT — WORKING 
CHAMBERS” ASSEMBLY 


Torque M,, indicated in the rotational hydraulic motor 
working chambers must be greater than torque М, loading 
the motor shaft (torque required by the driven machine 
(device)) because of the necessity of balancing also the 
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torque M,,,, of mechanical losses in the „shaft — working 
chambers? assembly. The assembly connects the shaft 
with working chambers, forms the working chambers and 
changes their capacity. Therefore, the torque М, indicated 
in the motor working chambers is a sum of torque M,, loading 
the shaft and of torque M „„ of mechanical losses: 

My, = Mu + Mum (1) 

Torque M,,,, of mechanical losses in a rotational hydraulic 
motor with variable geometrical capacity Чу Рег опе shaft 
revolution is, with maximum value of ae Le. with ae 
q, (with coefficient b, = = due Im =1), equal to the torque 
of mechanical losses in the motor working with constant 
theoretical capacity q,, per one shaft revolution. The theoretical 
and mathematical models describing the torque M,, of 
mechanical losses in the „shaft — working chambers” assembly 
of a motor with variable capacity q,,,, per one shaft revolution 
(with variable Б, coefficient) may be ‘described with reference 
to models of the torque Mum of mechanical losses in the 
assembly of a motor with constant capacity q,,, per one shaft 
revolution (with Б, = 1). 

Torque M,,. of mechanical losses in a rotational 
hydraulic motor is mainly an effect of friction forces 
between elements of the „shaft — working chambers” 
assembly, dependent, among others, on the torque M,, 
loading the shaft. 

Friction forces between the elements of the ,,shaft 
— working chambers” assembly are, to some extent, an 
effect of loading those elements by inertia forces from the 
rotational and reciprocating motion, dependent on the shaft 
rotational speed n,, and on the motor capacity dmg per one 
shaft revolution (b,, coefficient). 


There are also friction forces between the ,,shaft 
— working chambers” assembly elements and the working 
fluid, dependent on the fluid viscosity v and on the shaft 
rotational speed n,, and also on the motor capacity Amg Рег 
опе shaft revolution (b,, coefficient). The impact of working 
fluid viscosity v on the friction forces between elements of the 
„shaft — working chambers” assembly and the working fluid 
is visible mainly in the piston hydraulic motors with casing 
filled with fluid. 

Torque M,, loading the motor shaft and the shaft rotational 
speed п, required by the driven machine (device) change in the 
(0<@y< Oy, 0€ Mm < Му мах) hydrostatic drive operating 
range. The kinematic viscosity v of working fluid (hydraulic oil, 
oil — water emulsion) changes in the v... Xv X v, range. 

M n,, and v are parameters independent of the motor and 
of the losses in that motor. In models applied to motor with 
variable capacity dus bu du per one shaft revolution, the 


change of е, (b,,) is assumed in ће 0 < Amey (0 S б< 1) 


миш 1S Of the 0.2 + 0.3 order. 
of mechanical losses in a hydraulic 


range, although in fact b 
Torque M 


МиМкрлу Буру 
motor operating with torque M,, and speed n,, required by the 
driven machine (device) (the motor having capacity д, per one 
shaft revolution (b,, coefficient) and fed with working fluid of 


viscosity v) can be described as a sum of torque My ining on ce 


and increase AM, tap yy” 


> torque М имуу=бллу љу 
mechanical losses occurs in an unloaded motor (when the 
torque M, required of the motor is equal to zero — М, = 0) and 


the increase AM,, ofthe torque of mechanical losses is 
шМурлу Буру 


an effect of loading the ,,shaft — working chambers" assembly 
elements by the increasing shaft loading torque M,,: 


AM 


= + 
Mm|Mypnypbypv Mm|My =l ny bv (2) 


T АМ, мм, уву 
In constructing the theoretical апа mathematical models 


describing the torque AM,, of mechanical losses in 
Муру 


a rotational hydraulic motor, an assumption is made that the 


increase AM,, of the torque of the mechanical 
m|Myp%yp?yp¥ 


losses in the „shaft — working chambers” assembly, as an 

effect of the increasing required torque M,, loading the 

motor shaft, is practically independent of: 

— required shaft rotational speed n, 

— value of the b, = q,,,,/4),, coefficient of capacity per one 
shaft revolution, 

— working fluid viscosity v. 


An assumption was also made in the proposed models, that 


increase AM, of the torque of the mechanical 
m|Myp" yy 


losses is determined at the speed п, = п, i.e. equal to the 
motor shaft theoretical speed, at the coefficient Ь, = 1 
(with Amg 7 ймы) and at v = v, i.e. working fluid reference 
viscosity v : 


AM 


Mm|Myjp ny bay 


= AM 


= fM = 3) 


MmMy.mypbyaby, M 


In constructing the theoretical and mathematical models 


describing the torque M,, of mechanical losses in 
тМурпур Буру 


the ,,shaft — working chambers” assembly it was also assumed 
that: 


— the impact of required shaft rotational speed n,, and 
value of the coefficient Б, = q,.. /q,, on the load of the 
assembly elements with inertia forces, 

— the impact of working fluid viscosity у, п, апар, on the 
friction forces between elements and working fluid, 


in consequence, the impact of n,, р, and v on the torque 
Му 0f mechanical losses in the motor сап be determined 
with the shaft loading torque equal to zero (M,, — 0), 
because the inertia forces of the assembly elements and 
friction forces between the elements and working fluid are 
independent of the torque M,, loading the motor shaft: 
Mm|My -nM PMY = (пу, bs v) (4) 
The above mentioned assumptions allow to describe 


the torque М ым a p у Of mecha-nical losses in the „shaft 
M’ M’ M’ 


— working chambers” assembly by a theoretical model in the 
form: 


МһМулуруу Е 


+АМ G) 


Mm|My.nyreoby 7 1 Уң 


МиМуү=0,пуу,Буу,у 


In a hydraulic motor with theoretical (constant) capacity 
Gy, per one shaft revolution (b,, = 1), operating with theoretical 
(constant) shaft rotational speed n,, and the working fluid 
reference (constant) viscosity у, the theoretical model 
describing the torque of mechanical losses in the assembly 
takes the form (Fig.1): 


Mm| My nob "n 


- AM (6) 


Mm|My.nyeoby 7 1 Vn 


Mm|My 4-0.) p5)4= 1 Ma 


The impact of inertia forces of the „shaft — working 
chambers" assembly elements performing the rotational 


and reciprocating motion on the torque M, on , „v Of 

ммм 
mechanical losses іп an unloaded motor тау be presented 
as a function of the motor shaft rotational speed n,, and as 
a function of geometrical capacity Ф, (b,, coefficient) per 
one shaft revolution. 

The impact of the friction forces between the ,,shaft 
— working chambers” assembly elements and the working 


fluid on the torque M ymm -on._p,_y Of the mechanical losses 
M "Www 


in an unloaded motor may be presented as a function of the 
working fluid viscosity v and as a function of the motor shaft 
rotational speed п, and of the capacity а, (coefficient b, ) 
per one shaft revolution. 

The proposed theoretical models describing the torque 


M ymm аъ. , Of mechanical losses in the ,,shaft — working 
мем M 

chambers” assembly of an unloaded hydraulic motor 
(with torque M,, = 0), with changing shaft speed n,,, with 
theoretical (constant) capacity q,,, (b,, = 1) or geometrical 
(variable) capacity q,,,, (b,,) of the motor per one shaft 
revolution and with changing working fluid viscosity v, 
take the form: 

e ina hydraulic motor with theoretical (constant) capacity 


du (bu = 1) per one shaft revolution (Fig. 2): 


Mm|My-0,nyby-Lv ~~ Mm[My-0,n bel v, 


Va a (7) 
+ AM 


= (M iM sob М а) a 
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where: 


AM 


= + 
Mm|My=0,ny,byq=LVn Mm|My70,n y, by lv, 


-M 


Мъјм=0,пм=0м=1 7, (8) 
n 

5 -M oa ср у у 
Mm|My,, =0, ny, =п m -bm =l; Vn Mm|My =0,n4 =0,bq =1,у„7 үү 

Mt 


=(M 


* ina hydraulic motor with geometrical (variable) capacity 
“меу (dy, = Ом Чм per one shaft revolution (Fig. 5): 


aym 


M 


Mm|Myu=0.nybmY V Mm|My=0,n m bm:Vn 


v 


n 


=M +AM 


Mm|My-0,n 70,70, v, Mm|My=0,n е9. ти (9) 


= (М + AM 


Мт|Мм=0,п y, bm m ven 
n 


Mm|My-0,n 70b Lv, 


where: 


AM 


= + 
Mm|My=0,ny,.byVn Mm|My-0,n y, by. v 


Mm|My=0.n MERE (10) 


Mm|My-0,n y, by. v; B MM Oba Ly, 

n 
-M b 
Mm|M,,=0,n y=0,bųy=1,v M 
[My M M Jn, 


-(M 


Mm|My=0, nMn m -by=l Vn 


It is assumed in expressions (9) and (10) that the torque 


мым, 0.0, 0,4, of mechanical losses in the „shaft — working 


chambers" assembly of an unloaded hydraulic motor (M,, = 0) 
with geometrical (variable) capacity д. (Amey = Dy, Чу) per one 
shaft revolution and with zero shaft rotational speed (n,, = 0) is 
practically independent of the coefficient b,, of capacity per one 
shaft revolution and can be determined for b,, = 1. Therefore 
the following simplifying assumption is accepted: 


Mm|My4=0,n4=0,by4=0,v, 77 


= М 


Mm|My4=0,n)4=0,by4=1 Yn 


Mm|My-0,nyí70,byrvn = 


In models describing the torque of mechanical losses in 
a rotational hydraulic motor, used for description of losses 
and energy efficiency of the motor and of the hydrostatic 
transmission system in the (0 € Om < Фу, 0 < Mm < Mya) 
operating field, the increase of torque of mechanical losses 
occurring in fact at the motor speed n,, approaching the zero 
value (nų = 0) is not taken into account. That increase occurs 
below the shaft critical speed п. The motor rotational speed 
instability ón,, increases below the critical speed п, and 
in effect the torque M,, of mechanical losses in the ,,shaft 
— working chambers" assembly increases. Assessment of 


the value of torque Mani, озо, of mechanical losses 


in a motor with zero rotational speed п, (n,, = 0) is done by 
approximation of the function M = f (ny) at 
n, = 0. 

Exponent a, in expressions (7) and (9) determines the 
impact of the ratio v/v, of viscosity v to the reference viscosity 
v, = 35mm/’s" of the working fluid on the value of torque 
of mechanical losses. The impact occurs mainly in a piston 
displacement machines with fluid filling the casing (in pumps 
and hydraulic motors). 

The proposed theoretical model of the increase of 


torque AM,, of mecha-nical losses in a rotational 
mMm nM PMY 


Mm|My -0,ny by rn 


hydraulic motor, the increase resulting from loading the 
motor shaft with torque M,,, takes the same form in the 
case of a motor with theoretical (constant) capacity q,,, per 
one shaft revolution (b, = 1) and in the case of a motor with 
geometrical (variable) capacity Amg per one shaft revolution 
(dy, = by qud: : : " : 

e ina hydraulic motor with theoretical (constant) capacity 


du (b, = 1) per one shaft revolution (Fig. 1 and 3): 
AM (11) 


My 


а а Msi editado M 
Mt 


Mm[My, ny by =l, v 7 
=(M 


* ina hydraulic motor with geometrical (variable) capacity 
Sus bu auy per one shaft revolution (Fig. 4 and 5): 


AM (12) 
My 


elie oy euro M 
Mt 


Mm|My,ny,by,v 


-(M 


In effect, the proposed theoretical models describing the torque М, of mechanical losses in a hydraulic rotational 


motor take the forms: 


e ina hydraulic motor with theoretical (constant) capacity qu, (b, = 1) per one shaft revolution (Fig. 3): 


M 


Mm|Myq.ty.by=Lv 


i Mm|My 2 My; -n4 9n yi by Bic 


e ina hydraulic motor with geometrical (variable) capacity Qus (dy, = bu dug per one shaft revolution (Fig. 4 — with n 


Fig. 6 — with v... Vis Vmax): 


Мм, ambay  [УМи|Мм=0,ам=б, bL, Г (м 


+( МщМу=Муү›лм=пуү„бм=1,Уд - 
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Муму =0,пм=0,5м pu (М.м, =0,пм=пмиВм=1,% Муму =0,пм=0,5м a) 


Mm[My-0.n ny Ымм - 


aym 


Ny V 
Am Va 
| ) Mu (13) 
Mm|My =0,0y 2ny by =L; v, 
| My, 
Mt? MS 
aym 
уам ъ |У 
Mm|My=0,ny4=0,by=Lv, M 
Dt n 
(14) 


My 


Mm|My=0, ny =" bual) 
Mm 


My | Пы» Ям (b, -1 ) Yn 


| ИШ 


: | 
0 _ аР 
Ma OTT à 
Fig. 1. Torque M, of mechanical losses in the „shaft — working chambers” assembly of a rotational hydraulic motor with constant capacity q yy, 


Мт|М п, byl v 


n 
per one shaft revolution (b,, = I), at the shaft theoretical rotational speed n, 


and at the working fluid reference viscosity v, , as a function of the motor shaft 


torque M,,— graphical interpretation of the theoretical model (6) 


м=р — Мы=®,гм=0, №, м qu (6,21) 


1 | 
0 Vmin 4 Vmax ЕД 
V, Y, ГА 


Fig. 2. Torque М, 
Mm|M On bp bv 


of mechanical losses in the „shaft — working chambers” assembly of a rotational hydraulic motor with constant capacity q,, 


per one shaft revolution (b,, =1), at shaft torque М, = 0, as a function of the ratio v/v, of viscosity v to the working fluid reference viscosity v, — graphical 


interpretation of the theoretical models (7) and (8); motor shaft speeds: n,,= 0, п, n 


w The impact of working fluid viscosity v on friction forces between the 


„Shaft — working chambers " assembly elements and the working fluid occurs mainly in the piston hydraulic motors with casing filled with the working fluid. 


MATHEMATICAL MODELS OF THE 
TORQUE OF MECHANICAL LOSSES 


In mathematical models describing the torque М, оѓ 
mechanicallosses in a hydraulic rotational motor coefficients 
k, of losses are used relating (comparing) the components 
describing the torque М, of losses in theoretical models 
to the following reference values: 


— theoretical torque My = 1м т: of a hydraulic motor with 


theoretical (constant) capacity q,,, per one shaft resolution, 
determined at the decrease of Ap,, of pressure in the motor 


equal to the system nominal pressure р (Ap,, = p,) and with 
assumption that there are no pressure and mechanical losses 
in the motor, 

theoretical rotational speed п, of a hydraulic motor with 
constant capacity q,,, per one shaft revolution resulting 
from the motor capacity Q,, equal to the theoretical pump 
pt 

dat) 

theoretical capacity О of the hydraulic motor driving 
pump — a product of the theoretical capacity qp, per one 
shaft revolution ofthe constant capacity pump and the shaft 
speed n, of an unloaded pump (О, = q,, n,,). 


О, [ny = 
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ny=9, Пы м (by z1 ) Упіт, Fn, Ymax 


M 
Mim ч " = Mu pe m ^ zB _. 
ui SAT EORR 
\ ГА 


My © М 

пу = А f 
by=1 AMum ьи 
ГА Pa 


0 


Fig. 3. Torque Mma pny bgt of m echanical losses in the „shaft — working chambers” assembly of a rotational hydraulic motor with constant capacity q 


Ма= Pn М 


M 


per one shaft revolution (b, =1), as a function of the motor shaft torque M,,— graphical interpretation of the theoretical model (13); motor shaft speeds: 
theoretical speed п, n,,— 0; working fluid viscosity у, v, апау, „. The impact of working fluid viscosity v on friction forces between the „shaft — working 
chambers” assembly elements and the working fluid occurs mainly in the piston hydraulic motors with casing filled with the working fluid. 


My і пм 0 (by=0), Aug (by), Qus Qu (by=1 Jaa 
Mym ри 
» My M, =0 My 


Yn 


- Dy / n 
Mu йу = Mim pet AMum prc; 


Pa 


My We My 


Fig. 4. Torque М, мм. n © mechanical losses in the „shaft — working chambers” assembly of a rotational hydraulic motor with variable capacity 


MM Mn 


b 


due 


d, Per one shaft revolution, at theoretical shaft rotatio-nal speed n,, and at working fluid reference viscosity у, as a function of the motor shaft 


torque M,,— graphical interpretation of the theoretical model (14); capacity qw per one shaft revolution (coefficient b „of the change of capacity per опе 
shaft revolution): dus, = 0 (b,, = 0), d ue (b), Qus, 7 dy, (= 1) 


Theoretical and mathematical models describe the torque 
M um of mechanical losses in a rotational hydraulic motor with 
theoretical (constant) capacity q,, per one shaft revolution or 


with geometrical (variable) capacity usi bu du per one shaft 
revolution: 
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= is a theoretical capacity per one shaft 
9м Cup i-o; n растур 


resolution of a constant ca-pacity motor, determined at 
Ар = 0, p, = 9 and v,, which is equal to the theoretical 
active volume of the working chambers during one shaft 
revolution, 


n 
Mum by і 
y 
Mi 
Mum [bez 
Yn 
M M 
Mm 
by 
| T - pe 
0 І Vmin 1 Ymax » 
ГА Yn ГА 
Fig. 5. Torque М, of mechanical losses in the ,,shaft — working chambers” assembly of a rotational hydraulic motor with variable capacity 


Mm M 0.0 yp yp? 
du by Iq per one shaft resolution, at the shaft torque М, = 0 as a function of the ratio v/v, of viscosity v to the working fluid reference viscosity ү, 
= graphical interpretation of the theoretical models (9) and (10); motor shaft speeds: theoretical speed n p п, = 0; capacity q,,,, per one shaft revolution 
(coefficient b, of the change of capacity per one shaft revolution)): q,,,, = 0 (by, = 0), due = (b, i), Que, = dus (b,, = 1). The impact of working fluid viscosity v 
on friction forces between the „shaft — working chambers" assembly elements and the working fluid occurs mainly in the piston hydraulic motors with casing 
filled with the working fluid. 


"d n,70, ny, qu, 70 (buz0).quy, (bu) Iug Qu (bu 1), minh, Imax 

Mum бы 
м 

y 


| К 
Ям Pn 
0 Mu 2T M, 
Fig. 6. Torque M. _of mechanical losses in the „shaft — working chambers” assembly of a rotational hydraulic motor with variable capacity 


Mni M pn yp? y” 
dus, = b a4, per one shaft revolution as a function of the motor shaft torque M,,— graphical interpretation of the theoretical model (14); motor shaft speeds: 
theoretical speed т, п, = 0; capacity 4м, Рег one shaft revolution (coefficient b,, of the change of capacity per one shaft revolution): q ene 0 (b, 0), 
dus, (b, q, ui^ dus (b,, = 1). The impact of working fluid viscosity у on friction forces between the „shaft — working chambers” assembly elements and the 
working fluid occurs mainly in the piston hydraulic motors with casing filled with the working fluid 


— dug 0,915 а Geometrical capacity per one shaft revolution of the variation of capacity per one motor shaft revolution 
ofa variable capacity motor, determined at Ар, = 0, p,,,7 0 changes in the 0 <b,, < 1 range. 
and v „ which is equal to the geometrical active volume 
of the working chambers during one shaft revolution. The proposed mathematical models describing the torque 
In developing the models it is assumed that capacity Mm of mechanical losses in the „shaft — working chambers" 


Ivey per one hydraulic motor shaft revolution changes assembly, related to the theoretical models of the torque of 
in the 0 € 9м, < Чы range and coefficient b,, = 9а mechanical losses, take the form: 
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e inahydraulic motor with theoretical (constant) capacity qw 
(b,, = 1) per one shaft revolution [referring to theoretical 
model (13)]: 


Ny 


Mm|My,ny.by=lv— аманет Ми | + 
Mt n (15) 
n q р aym 
= M Mtt n 
+k,,.My = koi) tk, 2H pm +k, My 
Dy M? 
where: 
Ё _ M sink: оладат, e Муму eiiis HAW, 
FAS, 7 
Мм ЯмРь 
20 (16) 
k ‚Му ы йен Бы “mnie leu = 
7127 =Z 
Mm 
(17) 
M EE И и TR 
|. —— Mm|My,, -0.ny =n m by lv Mm|My, 20, ny -0,by =l; v. 
Ом‹йһ 
211 
k= Мар ме е МЫ ав „рм, 
727 M a C CECI 
Mt 
(18) 
M " "— € E 
= Mm|My, =Мм -NM 9 ny by =l; Vn Mm|My 20,ny 2n y; by zl Vn 
ЯмРһ 
2П 


• inahydraulic motor with geometrical (variable) capacity 
“меу (Ayer = Ом Au) per one shaft revolution [referring to 
theoretical model (14)]: 


vm 


n у 
= M 
Мула г bu My, p 
Mt n (19) 
aym 
n ймы | V 
= M Mt 
+k, .My= k; rt Kio bum emm mes tk; M, 
ia ӘП iv, 


where: coefficient k,,, is described by expression (16), 
coefficient k,,, — by expression (17), coefficient k,, — by 


expression (18). 
CONCLUSIONS 


]. Theoretical and mathematical models have been developed 
of the torque M,, of mechanical losses in the ,,shaft 
— working chambers" assembly of a rotational hydraulic 
motor with constant q,,, (V,,,) and with variable ds (V 
capacity per one motor shaft revolution. 

The models present dependence of the torque Mym of 
mechanical losses in the assembly on the torque M,, 
loading the shaft and on the shaft speed n,, [changing in the 
(0€ My Му 0€ ©m< мах) motor (and the hydrostatic 
transmission system) operating range] and also on the 
working fluid viscosity v changing in the v in S V < V range. 
Torque M,,, shaft speed n,, and working fluid viscosity v are 
independent of the motor and of the losses in it. 

The models present also the dependence of torque Myn 
on the capacity Amey per one shaft revolution (coefficient 


М, " 
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bu = qu, /q,, of the capacity per one shaft revolution) in 
a variable capacity motor. In the models, the change of Qus 
(б) is assumed in the 0 < Ss £ Im (0 < b,, < 1) range, 
although during the motor operation b,, 15 of the order 
of 0.2 + 0.3. 


2. Mathematical models of the torque М, of mechanical 
losses are based on the defined coefficients k, of energy 
losses, relating the torque of mechanical losses to the 
reference values: 

— theoretical torque M,, of the motor with theoretical 
(constant) capacity qu, per one shaft revolution 
determined at the system nominal pressure p, 

— theoretical rotational speed n,, of the motor with 
theoretical (constant) capacity q, , per one shaft revolution 
resulting from the pump theoretical capacity Q,. 


3. Mathematical models of the torque М,, of mechanical 
losses in the „shaft — working chambers” assembly should 
correspond with models of volumetric losses in the motor 
working chambers and with models of pressure losses in 
the motor channels. 
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ABSTRACT 


This paper presents results of research on spectral structure of underwater noise acoustic field radiated 

into water by selected ships moving in shallow waters. Underwater acoustic field of ships in motion is 

associated with acoustic activity of ship mechanisms and equipment under work. Vibration energy radiated 

by the mechanisms and devices is transmitted by ship structural elements to surrounding water where it is 

propagated in the form of acoustic waves of a wide frequency band. In this paper are presented results of 

the research on propagation of energy of acoustic waves in the near fiel, obtained from acoustic pressure 
measurements by means of two sensors located close to each other. 


Keywords: energy, propagation, hydroacoustics 


INTRODUCTION 


Research on spectral structures of underwater noise 
generated to water by vessels of different classes have been 
used for many years in the detection, location and identification 
systems, and also in the homing guidance and initiation units 
of sea weapons systems. 

Underwater sound field ofthe vessels is connected with the 
activity of acoustic wave sources installed on ships, i.e. ship 
mechanisms and equipment (main engines, generators, gears, 
pumps, shaft lines, pipes, ducts, etc.) and hydrodynamic sources 
such as screw propeller and water flow around the hull. 

Theoretical estimation of the level of acoustic energy 
near the wave source and attempts to present the vertical and 
horizontal distribution of energy in the plane of trajectory 
are difficult and currently not very effective. This follows 
from the fact that the theoretical model of extended source, 
which is a superposition of several different types of sources 
of acoustic waves, namely those due to water flow (propeller 
flow, flow around the hull), ship hull vibration and cavitation, 
makes that their full description is extremely complicated 
and difficult to apply in practice. A great difficulty is the 
necessity to take into account the influence of surface limiting 
water environment (free surface and bottom of sea). The 
next difficulty is the fact that the waves generated by a ship 
to water near the sources are non-stationary and nonlinear. 
It seems that many unknown variable parameters which 
shape acoustics of the near field of ships make application of 
adequate numerical models to solve this problem, impossible. 
For this reason, in practice, the measurements are performed 
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in natural conditions at various depths of ship operation and 
various set of points of their propulsion systems, and then 
results of the measurements are recorded and archived in 
digital form in specialty systems. 

To-be-archived sound pressure measurements should be 
recorded in the far field zone [1, 2] i.e. where the following 
condition is satisfied: 


d? » т)/2т or Ar « X (1) 


where: 


г  - the distance between source and receiver, 


Ar - distance between sensors, 

d - maximum linear dimension of the active part of 
transmitter, 

A  - acoustic wave length. 


Forthe assumed depth ofthe coastal water region, amounting 
from 10 to 70 m at which underwater noise measurements are 
made, the condition of the far field would be satisfied for the 
frequency f 31000 [Hz] in the trial area situated at the depth 
of 10 m. This limitation is also dependent on radiating surface 
dimension and sound propagation velocity in water. 

The sound pressure measurement performed on the 
measurement depths, does not allow to precisely characterize 
acoustic field distribution of vessel in motion. Characteristic 
components of the spectrum generated by the main and 
auxiliary mechanisms, propeller, shaft line of ships of different 
classes are contained in the frequency band up to about 150 
[Hz]. Sample distributions of hydroacoustic fields of selected 
vessels are shown in Fig. 1. 
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Fig. 1. Example acoustic spectrograms and power spectra of vessels, where: 
1 - ship with turbine propulsion system, 2 - ship with conventional propulsion system, 
a - the spectrograms made at the frequency band up to 500 [Hz] with resolution of 1/24 octave, 
b - the underwater noise spectrum made in the place where the most noisy part of the ship was situated directly above the acoustic sensor 


In the presented figures it is clearly visible that the 
characteristic spectral components, reflecting both the 
frequency and level of work of the main propulsion systems 
and auxiliary mechanisms, are located in the frequency band 
up to about 100 [Hz]. 

When taking into account the short characteristics of 
the acoustic field distribution of the above mentioned ships 
of different classes, it seems evident that that measurement 
ranges in the frequency domain should cover the whole range 
of length of acoustic waves propagating from the ships to 
surrounding water. 

Sound intensity measurement of underwater noise enables 
to carry out more accurate research in this water region. Such 
measurement based on a suitable spacing of acoustic sensors 
and precise phase matching of hydrophones and measuring 
setups, will make it possible to process of recorded data in the 
frequency band of interest. 

In order to determine spatial distribution of acoustic 
energy it is necessary to know vibration velocity distribution 
of acoustic wave sources. In the case of vessels to determine 
precisely vibration velocity distribution on hull shell plating 
Is very difficult because of a complex form of hull surface and 
limited number of measurement points. In association with the 
technical and financial difficulties, in the available literature [3, 
4] can be found a method for determining the acoustic wave 
velocity which makes it possible to determine a substitute 
vibration velocity distribution. Such distribution is determined 
experimentally by measuring the acoustic pressure and velocity. 
Measurement of acoustic pressure in an arbitrary point in space 
is easy to realize. For the measurement it is sufficient to process 
data from one sensor only, however the system consisted of at 
least two sensors is necessary to measure the acoustic velocity. 
The approximate component of acoustic velocity vector, 
directed along the vector r determined from acoustic pressure 
measurements with the use of two sensors placed close to each 
other, can be determined from the following relation: 


у= -—— |(рь -Padt 2 
r pAr B A ( ) 
where: 

P- р, - pressure difference at the points A and B, 

р - density of medium. 


The approximation can be applied to the far field in which 
the relation (1) is satisfied. During the acoustic velocity 
measurements the measured values of sound pressure are often 
converted from analogue to digital form. 


van- ip, @-p,@) O 
pAr 


where: 
v, - constant component of velocity which can be determined 
after calculation of velocity values for the entire wave 
period. 


The resulting acoustic velocity is used to determine the 
acoustic wave intensity which, for stationary signals, can be 
described by the dependence: 


I(r) = p(r.t)v(r.t) (0 
Taking into account ће dependence (3) we obtain: 
| —————— 
I^ zo Pa * Ps) |(р» =рл) (5) 


And, after an appropriate transformation the digital 


dependence (5) takes the form: 
1 

I= ——— - In JG, (k 6 

oA [G,, (K)] (6) 


where: 
0<k<N 
Im[G..(k)] - complex part of cross spectrum. 
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LIMITATIONS OF THE PRESENTED 
METHOD 


The method used to measure sound intensity with the use of 
two sensors located close to each other, has some limitations. 
The basic error of this method is connected with inaccuracy 
of measuring the derivative of hydroacoustic pressure. The 
error is observed in the higher frequency ranges. For a flat 
harmonic wave, the estimated value of related to the exact 
value of { the sound intensity I, can be determined from the 
dependence (7): 


1 _ sin(kAr) (7) 
1 kAr 
where: 
k - wave number. 
In the logarithmic form the dependence is as follows: 
L=10 "E (8) 
kAr 


For the spherical wave, the estimated value of Í related to 
the exact value of the sound intensity I, can be determined from 
the dependence (9): 


i_ sin(kar) d Е | ` D 


I kAr 4 


An essential error of the method in the low frequency ranges 
is a deviation resulting from the phase mismatch of measuring 
setups [2, 3]. For this reason, we should tend to exactly match 
the phases of all the measuring setups at which ọ — 0. If the 
phase mismatch occurs between the two setups for flat wave 
then the estimate of | - value related to the exact sound intensity 
I, can be determined by using the relation (10): 

I _ sin(kAr+@) (10) 
1 kAr 
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CALIBRATION OF MEASURING SETUPS 


In order to obtain authentic results of measurements it is 
necessary to perform accurate calibration of measuring setups 
to determine the exact phase difference of measuring setups. To 
determine the differences the method of submission of mutually 
perpendicular harmonic signals has been used. The curve 
obtained as a result of the composing of signals is dependent 
on frequency and phase of the signals as well as values of their 
amplitudes. Results of the composing of signals were published 
first by Lissajous, a French scientist, at the end of 19" century. 

The laboratory calibration tests of the measuring setup of 
the probe for underwater, measuring the sound intensity have 
been performed in a measuring tank. The block diagram of the 
measuring system including the tank, is shown in Fig. 2. 

In order to eliminate the waves reflected from the tank walls 
and water surface the system has been applied for calibrating 
the sensor by means of the impulse method [2], which makes 
it possible to eliminate unfavorable reflected waves by suitable 
selection of duration time of pulses of harmonic waves generated 
in appropriate time intervals. The duration times of sending 
pulses resulted from the basin’s geometry and were determined 
for the reflections from the tank walls from the dependence: 


Мт? +b? -r 


- for reflections from the bottom and surface: 
2 2 
vr —-h*-r 
t < —————— (12) 
Cc 
- for the reflections between hydrophones: 


(11) 


r 
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(13) 
where: 
T - 
с - 


duration time of a signal free from reflected signals, 
sound velocity in water. 


Fig. 2. Block diagram of the measuring system including the tank. Notation: 1 - sending hydrophone, 2 - receiving hydrophones 
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From the far field condition (1.1), by knowing the maximum 
linear dimension of the transmitter, d < 0.1 m, and the sound 
velocity in water, c = 1480 m/s, the minimum distance from the 
receiving hydrophones to the sending one, has been determined 
at the frequency for which the measurement system has been 
controlled. 


Tab. 1. The length of the near field transmitter for 


f [kHz] 5 6 7 8 9 10 
d^f/c [m] | 0.033 | 0.041 | 0.047 | 0.054 | 0.061 | 0.068 


The measurements were performed for the distance between 
receiving and sending hydrophones, r < 0.4 m, and the pulse 
duration which could not be longer than t < 0.4 s, taking into 
account the dependence (11), (12) and (13). 

The tested probe consisted oftwo 4032 RESON hydrophones 
spaced at the distance Ar = 38 mm, and having the sensitivity 
values: 

A - no. 3702080: 170.4 dB at 250 Hz ; 169,1 dB at 8 kHz, 
B - no. 3702082: 170.3 dB at 250 Hz ; 169 dB at 8 kHz. 


The difference in sensitivity of the sensors in the band of 
interest, up to 10 kHz, is negligible. The significant differences 
in the indicated levels, reaching from 3 to 10 dB, occurred in 
the frequency band from 20 to100 kHz. The 8100 Brüel & Kjar 
hydrophone which generated acoustic wave in the frequency 
ranges given in Tab. 1, served as a source of harmonic waves. 
The sinusoidal signal from the generator was transmitted 
to the 4440 Brüel & Kjar gate which controlled the length 
and pulse repetition frequency. The signal from the gate was 
properly amplified by 2713 Brüel & Kjar power amplifier 
and then transmitted in the form of short pulses of harmonic 
wave, to the surrounding water by the sending hydrophone. 
The acoustic wave produced by the hydrophone was received 
by two receiving hydrophones placed at equal distances from 
the sending sensor. Signals received by them were amplified 
by 2636 Brüel & Kjar amplifiers, and then recorded and 
processed in a analyzer which made it possible to determine 
mismatch of the tested phases of the system. The extensive 
results of the performed calibration of the probe are contained 
in the report from realization of this work [2] performed at 
the Radiolocation and Hydrolocation Department of Polish 
Naval Academy in Gdynia. The publication presents only 
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representative results of the probe calibration, recorded in the 
form of Lissajous curves. 

From the measurements performed for the probe in 
question it can be concluded that the phase mismatch of the 
measuring setups did not exceed 0.2? for the tested wave 
length. The level of phase mismatch was enlarged to 0.3? after 
taking into account the specific conditions of measurements 
conducted in the dynamic trial areas. The lower frequency 
range for different distances between sensors was calculated 
for this value (acc. Eq. 10). 
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Fig. 4. The values of the phase mismatch level of the lower frequency range 
for the distance between sensors Ar = 0.2, 0.4, 0.6, 0.8, 1, 1.25, 1.5 m, 
respectively 


The upper frequency range for the same distance was 
determined in the same way by taking into account the 
dependence (7). On the basis of those frequency ranges the 
frequency range at which the error does not exceed 1 [dB], 
was determined. 

On the basis ofthe performed research it can be concluded 
that, for the phase mismatch level of 0.3?, measurements in 
dynamic trial areas can be conducted in the frequency range 
from 10 [Hz] +2 [dB] to 1325 +1[dB] and at the spacing 
between hydrophones equal to 0.2 [m]. The described sound 
intensity measuring method makes performing the tests within 
the specified frequency band in the coastal zone, possible. 
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Fig. 3. The result of probes calibration matching 
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Fig. 5. The frequency ranges of the sound intensity measurement at which 
the error does not exceed 1 [dB] 


DYNAMIC MEASUREMENTS 


The acoustic noise measurements of submarines are 
repeated periodically during operation of each ship. During the 
measurements ships pass at least twice through the trial area 
with the set work parameters of the propulsion system. The 
determined ship parameters are reached at the distance of 300 m 
at least before the trial area and maintained over the distance 
of 600 m at least (300 m behind the buoys). Sound pressure 
measurements are made at a distance afore and astern the ship 
to make it possible to characterize underwater disturbances of 
the ship. The passage of ship through the trial area is very often 
presented in the form of sound pressure in function of time 
and frequency. The spectrogram of the performed recording 
is presented in Fig. 6. 
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Fig. 6. The acoustic field spectrogram 
of the ship moving with forward speed of 4 kn 


The spectrogram consists of 299 spectra recorded every 312 
[ms] with the resolution of 0.04167 octave in the frequency 
range from 3[Hz] to 2.818[kHz]. Broken lines mark the distance 
of the ship from the acoustic sensor. 

In the spectrogram two areas are clearly visible. The 
first area contains frequencies to about 100[Hz]. In this 
area the characteristic spectrum components resulting from 
working ship mechanisms, can be observed. The second area 
of the frequency range from 100[Hz] to 2.8[kHz] contains 
a continuous spectrum. The spectrum is connected with work of 
cavitating propeller, turbulent flow in pipelines, flow around the 
hull, airflow in ventilators, etc. Such analysis makes it possible 
to select, out of the spectrogram, a spectrum in an arbitrary area 
covered by the measurements. 

Fig. 7 presents the ship’s spectra taken from the places 
pointed in the spectrogram. The first of them was recorded when 
the most noisy part of the ship was situated about 50[m] before 
the acoustic sensor. The second - when the disturbance source 
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was just above the sensor, and the third one - corresponds with 
the distance of 50[m] behind the sensor. 

The two previously discussed areas are also visible in the 
spectrum in question. From them it can be also observed that 
before the trial area the ship generates mainly waves connected 
with work of mechanisms and ship equipment, and behind it 
the waves connected with work of the propeller are dominant. 
To more clearly present the figures, the signal recorded when 
the ship passed just above the acoustic sensor, was suppressed 
by 30[dB], and that corresponding with the distance of 50[m] 
behind the trial area - by 60[dB]. 
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Fig. 7. The spectra obtained at some distance before and behind the trial 
area, as well as during passing just above the acoustic sensor. 
Notation: 1 - the spectrum of the ship taken 50[m] behind the trial area, 

2 - the spectrum recorded when the most noisy part of the ship passes above the 
acoustic sensor, 3 - the spectrum of the ship taken 50[m] before the trial area 
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The characteristic components visible in the first area 
can be unambiguously assigned to working mechanisms and 
equipment of the ship. Examination of ship spectrum is one 
of the methods for identification of ship underwater noise. 
In order to carry out identification of the components such 
tests are usually performed in two stages. The first contains 
measurements carried out on anchored ship. The measurements 
consist in the measuring of vibration of main mechanisms 
and auxiliary equipment of the ship, accompanied with the 
simultaneous measuring of sound pressure in water depth. In 
the second stage the ship crosses dynamic trail areas, under 
various operational settings of its propulsion system. During the 
measurements both underwater noise and vibration of selected 
mechanisms and equipment, are recorded. On the basis of the 
measurements it is possible to characterize the components 
visible in the first presented field. To identify the components, 
an analysis with using constant band width filters for the 
tested frequency band, is often performed. A representative 
identification analysis of the spectrum of the ship in motion is 


shown in Fig. 8. 
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Fig. 8. The acoustic field spectrogram of the ship with forward speed 
of 4 kn, recorded in the frequency band up to 100[Hz] 


The spectrogram consists of 74 spectra recorded every 
1.333[s] with the resolution of 0.25[Hz] in the frequency band 
up to 100[Hz]. From the figure was selected the spectrum 
when the ship's power plant was situated just above the 
acoustic sensor (the place distinguished with black line at the 
spectrogram). 

The figure clearly shows discrete spikes coming from 
working main engines, shaft lines and propeller, and a single 
spike marked “Т” in red color, coming from an electric generating 
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Fig. 9. The ship acoustic power spectrum recorded at the distance marked “0” in Fig. 8 


set under operation. Each discrete spike was distinguished by 
a successive number, namely green numbers indicate spikes 
coming from propeller and shaft lines, and black numbers 
show spikes coming from working main engines. Additionally, 
each discrete spike was distinguished by broken line, namely 
blue lines stand for the frequencies connected with work of 
main engines, and black lines - the frequencies connected 
with rotation of shaft line and propeller. The identification 
results obtained by using the method in question have been 
published many times [8, 12, 16]. This paper presents only 
a single analysis performed for a given working regime of 
ship propulsion system. The authors’ archive contains results 
of identification of ships of different classes, made for various 
working regimes of their propulsion systems. 

Among the relationships which are most commonly 
elaborated on the basis of such measurements, can be numbered 
the above presented characteristics of the sound pressure level 
in function of frequency, obtained for an arbitrary phase of 
ship’s passing through the trial area, as well as the following: 
a) the sound pressure level in function of distance, 

b) the sound pressure level in function of ship forward speed 
and water depth, 

c) the relation between underwater noise and vibration of ship 
mechanisms. 


The sample characteristic relation of sound pressure level 
in function of distance between ship and sensor, and of ship 
forward speed, is shown in Fig. 10 and 11. The curves were 
determined on the basis of many repeated measurements 
performed in the trial area located in the same place, and 
during the measurements of the ship crossing the trial area 
under the same work regime of its propulsion system. Results 
of the measurements were then subjected to relevant statistical 
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Fig. 10. The characteristic relation of the sound pressure level in function of 
distance between ship and acoustic sensor. Measurements of the ships were 
performed at the water depth h — 30 [m]. Notation: 

1 - the measurements performed at the speed v — 14.5[kn], 

2 - the measurements performed at the speed у = 12.5[kn], 

3 - the measurements performed at the speed v = 8.5[kn], 

4 - the measurements performed at the speed v — 5.5[kn], 

5 - the measurements performed at the speed v — 2.5[kn] 


processing and approximating procedure. The approximation 
by means of polynomial of the second or third order is most 
commonly used to approximate the relation of pressure levels 
in function of ship distance, water depth and ship forward 
speed. 
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Fig. 11. The relation of the sound pressure level as a function of ship s 
forward speed. Notation: 
1 - the measurements performed at the water depth h = 10[т], 
2 - the measurements performed at the water depth h = 20[m], 
3 - the measurements performed at the water depth h = 30[m], 
4 - the measurements performed at the water depth h = 40[m], 
5 - the measurements performed at the water depth h — 50[m] 


Inthe paper possible ways of extending the measuring range 
due to a change in measurement methods, are mainly presented. 
To make performing the sound intensity measurements 
possible, the trial area was additionally equipped with a pair of 
sensors tested in laboratory, and positioned along the assigned 
path of ship through the area. In order to verify correctness 
of indications of the probe, the measurement results of noise 
propagating from the sensors of the probe were compared to 
each other, and the hydroacoustic pressure levels in the selected 
frequency bands were also compared. 

The equal values of the levels obtained from the sensor 
probe, in the considered frequency bands, and given in the point 
(L = 138 [dB], f= 930,6 [Hz], t = 23,10 [s]) distinguished by 
the cursor at the spectrograms, and the character of underwater 
noise distribution confirmed the possibility of conducting 
mutual analyses by the tested probe, which should enrich in 
consequence knowledge on underwater noise. 

One of the fundamental features which differed the sound 
intensity measurement from the sound pressure one 1s the 
possibility of determining the phase differences between 
active and passive part of the sound field [3], which allows 
to determine direction of propagation of acoustic waves in 
water. 

In view ofthe limitations imposed on the application ofthe 
measurement method in question, revealed by the equations (10) 
and (12), the frequency band for processing the measurement 
results was reduced from 6.9 [Hz] to 1.334 [KHz]. The accuracy 
of value of the level for the frequency of 6.9 [Hz] should not 
exceed 3 [dB], and for the upper frequency: 1 [dB]. 
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Fig. 12. The measurement results of hydroacoustic pressure recorded by using the probe 
for measuring the sound intensity of the ship moving with forward speed of 6.6 kn. 
Notation: 1, 2 - the spectrograms of the ship acoustic field obtained from the sensor probe at the frequency band from 6.9 [Hz] to 11.22 [kHz], 
a, b - the relation of the pressure levels in function of duration time of ship passage through the trial area, obtained for the following frequency bands: 
І - the pressure value L = 155 [dB] in the band from 6,9 [Hz] to 11.22 [kHz], 
II - the pressure value L = 154 [dB] in the band from 100 [Hz] to 11.22 [kHz], 
ШІ - the pressure value L = 151 [dB] in the band from 6.9 [Hz] to 100 [Hz]. 


From the pressure measurements presented in the 
spectrogram of Fig. 12 the sound intensity was determined in 


accordance with the equation (1.6). 
[s] (Time) ^ Calc. Intensity Spectrum() (Magnitude) [dB/1,00p W/m 


— 110 
102 


[Hz] 


Cursor Values Status Delta 
Y = 74.6 dB/1.00p Wim? 2006-04-20 11:34:33.073 Sum = 112dB-/1.00p Wim? 


Averaging time: 1s 
X = 53.86 Hz 
2= 24.105 


Fig. 13. The spectrogram of underwater noise intensity obtained in the 
frequency band from 6.9 [Hz] to 1.334 [kHz] 


When testing the sound intensity amplitude of passage of the 
ship with the course of 0° relative to the location of the axis of 
the sensor probe, it should be demonstrated that on the basis of 
measurements in real conditions it is possible to precisely locate 
the place of passage of the ship above the set of acoustic sensors. 
To determine the location is possible because in laboratory 
conditions the value of sound pressure level falls to “0” for 
perpendicular position of the sound source relative to the probe. 
The phenomenon is described in many publications [4, 5, 6], and 
confirmed in laboratory tests carried out on the probe. A local 
depreciation of sound intensity revealed on the spectrogram are 
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visible for about 24 [s] duration time of passage of the ship above 
the probe. For different wave lengths we can observe of course 
different instants in which such local reduction of the sound 
intensity level would occur. The sound pressure level spectrum in 
function of time for the previously selected frequency ranges was 
produced to accurately reveal the local minimum values. In the 
spectra shown in Fig. 14, are clearly visible the local minimum 
values of the sound intensity level, by which perpendicular 
direction of hydroacoustic wave propagation relative to the probe 
in selected areas, is determined. 
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Fig. 14. The spectra of sound intensity level as a function of duration time 
of passage of the ship above the probe, obtained for the selected frequency 
bands: 1 - the minimum value of sound intensity level L = 93.4 [dB], for the 
time t = 22.90 [s], in the frequency band from 6.9 [Hz] to 1.334 [kHz]; 

2 - the minimum value of sound intensity level L = 94,7 [dB], for the time 
t — 22.30 [s], in the frequency band from 100 [Hz] to 1.334 [kHz]; 

3 - the minimum value of sound intensity level L = 91.3 [dB], for the time 
t = 24.10 [s], in the frequency band from 6.9 [Hz] to 100 [Hz] 


The places of occurrence of local minimum values of sound 
intensity level determine a change of direction of energy flow 
to the probe. The change of direction can be showed by using 
the analysis of real part of complex spectra of sound intensity 
level. Results of the performed analysis is shown in Fig.15. 
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Fig. 15. The spectrum of real part of complex spectra of sound intensity 
level, performed in the frequency range from 6.9 [Hz] to 1.33 [kHz] 


On the presented diagram of the analysis of absolute 
value of sound intensity level in function of duration time 
of passage of the ship above the probe the discussed change 
in energy flow direction is clearly visible. Positive value of 
sound pressure can be observed up to 22.90 [s] of duration 
time of passage of the ship through the trial area, represented 
by the black curve in Fig. 15. The black curve shows negative 
value of sound intensity. And, the maximum value of sound 
intensity L = +118 [dB] re 10? [W/m2] ( distinguished by 
cursor ) is observed at t — 21.10 [s], and the minimum value 

= -114 - at t = 26.80 [s]. The spectra obtained from the 
spectrogram presented in Fig. 13, which were taken for an 
arbitrary approaching time of ship to the probe and during its 
departure from the trial area, can also show the reversal of 
sign of sound intensity level. 

From the spectra presented in Fig.16 and 17, the summary 
sound intensity levels equal to L = +113 [dB] for the 
upper spectrum and L = -103 [dB] for the lower spectrum, 
respectively, can be read. The reversal of sign of sound intensity 
level constitutes a confirmation that reversal of direction of 
energy flow to the sensor probe takes also place. 

Reversal of direction of energy flow, determined from the 
measurment of difference of pressure values obtained from two 
hydroacoustic sensors makes it possible to determine direction 
ofthe source relative to location of the probe, within the sector 
from -90? to +90°. Determination of the source direction by 
means of this method is rather inaccurate. 

Such measurement informs only whether the acoustic noise 
source is situated before the perpendicular plane located in 
half-distance between sensors, or behind it. 


A more accurate bearing of the probe equipped with a pair 
of sensors can be obtained by testing both the sign of sound 
intensity and phase of the signal between the active and 
passive part of the sound field. For this purpose the tested noise 
was so processed as to obtain the spectrogram of amplitude 
of underwater noise phases, from which it was possible to 
determine the phase [°] within the considered frequency 
ranges. 
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Fig. 16. Spectrum of sound intensity in function of frequency, obtained 
during approaching the ship to the probe. 
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Fig. 17. The spectrum of sound intensity in function of frequency, obtained 
during departing the ship from above the probe 


The spectrogram of phase amplitude of underwater noise 
(phase-assigned auto-spectrum) was calculated from the 
dependence: 

Oxy 
Pyy = 


14 
" (14) 


where: 
x = the active component of sound field, 
y 7 the reactive component of sound field. 


The spectra of underwater noise phases were determined 
in the selected frequency bands, Fig. 18.III. They are only 
useful in determining the sign of phase during performance 
of the measurement. The spectra were determined from the 
underwater noise phases, Fig. 18.1, for the whole duration time 
ofthe measurement. Examining the figures we can observe that 
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Fig. 18. I - the spectrogram of phase amplitude of underwater noise, П - the spectrum of phases of underwater noise, Ш - the spectrum of noise phases 
calculated for the selected frequency bands: 1 - that from 6.9 [Hz] to 1.334 [kHz], 2 - that from 100 [Hz] to 1.334 [kHz], 3 - that from 6.9 [Hz] to 100 [Hz] 


during approaching the trial area by the ship to, at t = 23 [s], the 
phase sign was negative, and behind the area - positive. 

By comparing the spectra of sound intensity and signal 
phases it can be observed that during approaching the probe 
by the ship the sound intensity level obtains positive value and 
the phase - negative. The obtained positive values of sound 
intensity do not require any comments. The phenomenon is 
precisely described in many literature sources [1, 2, 4, 5, 6]. 

The auxiliary Fig. 19 can be used to explain how to use sign 
value for a more accurate location of sound sources. 

C 


Fig. 19. Summation of vectors of waves of different amplitudes and initial 
phases. Notation: x, — the deflection of the active component of the sound 
field wave amplitude r , initial phases д, and angular velocity co, 

x, — the deflection of the passive component of the sound field wave 
amplitude r „ initial phases д, and angular velocity о, 

x —the deflection of wave resulted from summation 
of component vectors of active and passive sound field 


From Fig. 19 wecan easily determine the phase ofthe signal 
by using the dependence: 
ies BC  mnsinó, +r sind, 


= (15) 
OB r,cosd, +r, cosd, 
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From principles of trigonometry it is known that tangent 
function takes positive values in the 1st and 3rd quadrant, and 
negative — in the 2nd and 4th one. 

Based on the knowledge of sign of sound intensity level and 
signal phase, a sound source can be located in an appropriate 
quadrant determined by the coordinates of sensor probe. The 
pictorial scheme of the location of the sound source in the 
relevant quadrant determined by pair of acoustic sensors, is 
shown in Fig. 20. 


course [0°] 
course [180°] 


Fig. 20. The location of ship in the relevant quadrant determined by the 
sensor probe coordinates, during ship 5 passage trough trial area. 
Notation: I - the area of positive values of sound intensity level, 

II - the area of negative values of sound intensity level, 

a - the quadrant of positive values of the sound field phase, 

b - the quadrant of negative values of the sound field phase, 

1, 2 - the hydroacoustic sensors 


After the measurements carried out for the course of 
0 [?] relative to the location of probe acoustic sensors, the 
measurements for the course of 180 [?] were performed. During 
the measurements the ship did not changed the area covered 
by the tests, against the location of the sensors axis, as shown 
in Fig. 18. From the tests the spectrogram of underwater noise 
intensity and spectrum of sound intensity in function of duration 
time of ship passage through the trial area, were achieved. In 
the spectrum it is clearly visible that when the ship approaches 
the probe (up to t = 43 [s]) sound pressure levels take negative 
values, and positive ones during the ship's departure from it. 
From comparison of the results of the measurements taken for 


the course of 0 [°] and 180 [°], the difference in signs of the 
sound intensity level and phase during approaching the probe 
by the ship and departing from it, is clearly visible. 
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Fig. 21. The spectrogram of underwater noise intensity recorded in the 
frequency band from 6,9 [Hz] to 1.334 [kHz] 
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Fig. 22. The analysis of the real part of complex spectrum of sound 
intensity, performed in the frequency band from 6.9 [Hz] to 1.33 [kHz] 


On determination of the sound intensity, the determining 
of signal phases whose results are presented in Fig. 23 and 24, 
has icta started. 
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Fig. 23. The phase amplitude sae sec ат of underwater noise, recorded in 
the frequency band from 6.9 [Hz] to 1.334 [kHz] 
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Fig. 24. The spectra of phase amplitude signals of underwater noise, 
recorded in the frequency band from 6.9 [Hz] to 1.334 [kHz] 


Comparing the results of phase measurements carried 
out for the tested courses, one obtained also the difference in 
signs, which confirms that to determine ship's course relative 
to sensors is possible. 


CONCLUSIONS 


The results of the performed research have demonstrated 
that by testing the sign of the sound intensity level and phase 
of noise signals it is possible to locate sound source just in 
a relevant quadrant determined by the sensor probe coordinates. 
To achieve a better accuracy is possible by increasing the 
number of pairs of sensors and spacing them properly. The 
use of three sensors located in an equilateral triangle will 
make it possible to reduce the sector of sound source location 
to 30 [?]. 

As demonstrated in the paper, the probe applied to measuring 
the sound intensity will make it possible to supplement the so 
far collected results of underwater noise investigations carried 
out in shallow waters by knowledge on direction of propagation 
of hydroacoustic waves. 

On the basis of known values of sound pressure and intensity 
it is possible to elaborate the following characteristics: 

a) the noise spectral structures of vessels, 

b) the results of vibration of mechanisms and equipment 
installed on ships, which greatly affect spectral structure 
ofship acoustic field (main engines, generators, shaft lines, 
screw propeller), 

с) the results of research on the relationship between ship 
underwater noise and vibration of ship mechanisms, 

d) initial assessment of technical state of ship machinery, 

e) estimation of the energy transmission coefficient, 

f) the measurement results of sound pressure level changes 
in function of ship's forward speed, 

g) the results of sound pressure level changes in function of 
distance, 

h) presentation ofthe acoustic wave propagation model made 
for water depth values selected along ship's trajectory, 

1) determination of direction of propagation of hydroacoustic 
waves. 


This paper has an experimental character and constitutes 
a continuation of many other investigations in the subject- 
matter area [7 +18]. 

In the future a series of extensive studies should be 
conducted, on the basis of which the maximum distance from 
which to identify and locate particular vessels is still possible, 
could be determined. The investigations should be performed 
in the coastal zone, at different sea states and different levels 
of hydroacoustic background noise, as well as in the trial areas 
located in the Gdansk Bay, at different water depths. 


BIBLIOGRAPHY 


1. Kozaczka E., Grelowska G., Bittner P., Baranowska A., 
Milanowski W., Dobrzaniecki J.: Spatial distribution of 
underwater noises radiated by ship to hemisphere and its 
standardization Gdynia 1997 

2. Kozaczka E., Grelowska G., Gloza I., Dobrzaniecki J.: 

Determination of phase displacement of hydrophone. Control 
measurements in a test tank and at sea, Hydroacoustics, Volume 
3, Gdynia 2000 

. Bulletin of Brüel &Kjar, 1973 

. Elliots S.J.: Errors in Acoustic Intensity Measurements. Journal 

of Sound and Vibration 1982 

5. Thompson J.K.& Tree D. R.: Finite difference approximation 
errors in acoustic intensity measurements. Journal of Sound and 
Vibration 1981 


+ о 


POLISH MARITIME RESEARCH, No 3/2010 35 


6. Jacobson F.: Measurement of sound intensity. Acoustics 
Laboratory of the Danish Technical University. 

7. Kozaczka E., Grelowska G., Bittner P., Baranowska A., Kicinski 
W., Milanowski W.: Spatial distribution of underwater noise 
radiated by ship to hemisphere. STAGE П, Proceedings of 9th 
Symposium on Hydroacoustics, Gdynia 1992 

8. Baranowska A., Gloza I.: Identification of underwater 
disturbance sources with using the coherence function. 
Proceedings of 48th Open Seminar on Acoustics, Wroclaw 2001 

9. Dobrzaniecki., Gloza I., Domagalski J.: Preliminary research 
on monitoring the own ship noise, Proceedings of 50th Open 
Seminar on Acoustics, Szczyrk — Gliwice 2003. 

10.Gloza I.: The investigation of changes of ships underwater noise 
in shallow sea, Doctoral thesis, Polish Naval Academy, Gdynia 
1994 

11.Gloza I.: Experimental research of noise generated by ship 5 
engine, Proceedings of 6th Symposium on Hydroacoustics, 
Gdynia 1989 

12.Gloza I.: Assessment of technical state of ships propulsion 
system on the basis of measurement of hydroacoustic noise. 
Proceedings of 10th Symposium on Hydroacoustics, Gdynia 
1993 

13.Gloza I., Malinowski S.: /dentification of ships underwater noise 
sources in the coastal region. Hydroacoustics, Vol. 5/6, 2003. 

14.Gloza I., Domagalski J.: The investigation of propagation of 
acoustic waves generated by moving ship. Hydroacoustics, Vol. 
5/6, 2003. 

15.Gloza I., Domagalski J., Malinowski S.: The character of 
underwater noise radiated by small vessels. Hydroacoustics, 
Vol. 4, 2001. 


36 POLISH MARITIME RESEARCH, No 3/2010 


16.Gloza I., Domagalski J., Malinowski S.: /dentification of ships 
underwater noise sources in near field. Proceedings of 49th 
Open Seminar on Acoustics, Warszawa - Stare Jablonki 2002 

17.Gloza I., Domagalski J., Dobrzaniecki J.: The structure of 
acoustics field of vessel s underwater noise in near field. 
Proceedings of 48th Open Seminar on Acoustics, Wroclaw 2001 

18.Grelowska G., Bittner P.: The processing of underwater signals 
produced by ship with the possibility of increasing the frequency 
resolution. Scientific Bulletin of Polish Naval Academy, vol. 
XXXIV, 4, 1994. 


CONTACT WITH THE AUTHORS 


Eugeniusz Kozaczka, Prof., 
Faculty of Ocean Engineering 
and Ship Technology 
Gdansk University of Technology 
Narutowicza 11/12 
80-233 Gdansk, POLAND 
fax: (058) 347-21-81, 
e-mail: kozaczka@pg.gda.pl 


Jacek Domagalski, Ph. D. 
Ignacy Gloza, Ph. D. 

Faculty of Navigation and Naval Waeapons, 
Polish Naval Academy 
Smidowicza 69 
81-103 Gdynia POLAND 
e-mail: i.gloza@amw.gdynia.pl 


POLISH MARITIME RESEARCH 3(66) 2010 Vol 17; pp. 37-44 
10.2478/v10012-010-0026-z 


Mathematical model 
of radial passive magnetic bearing 


Leszek Matuszewski, Ph.D. 

Gdansk University of Technology 
Krzysztof Falkowski, Ph.D. 

Military University of Technology, Warsaw 


ABSTRACT 


In the article a mathematical model of radial passive magnetic bearings applicable to ocean engineering 
units has been presented. The application of the bearings in ship thrusters should increase durability of 
propulsion systems and give lower maintenance costs. The rotor of thruster s electric motor is suspended 
in magnetic field generated by radial passive magnetic bearings. However the maintaining of axial 
direction of the rotor must be controlled with electromagnets equipped with a high-dynamic controller. 
The risk of application of the magnetic bearings results from very low stiffness of the unit in comparison 
with rolling bearings. Also construction of the bearing should be different due to gyroscope effect and high 
forces generated during ship manoeuvring. The physical model performs correctly if electromagnets are 
controlled properly; and, technological problem with sealing system seems more significant than power 
supply to electromagnets winding. The equations presented in the article are necessary to build algorithms 
of a digital controller intended for on-line controlling magnetic bearing in axial direction. 


Keywords: propeller, ring motor, magnetic bearing 


INTRODUCTION 


Possible application of magnetic bearings to ocean 
engineering units has been considered and tested for several 
years. Their low resistance to rolling and very long life time 
is of a great advantage. However because of a small stiffness 
they require introducing constructional changes into existing 
devices and new designed ones to be fitted with the so-called 
floating rotor. The bearing system presented in this paper is 
superior, as regards its durability, to traditional mechanical 
bearings both sliding and rolling ones. The most promising 
example of application of such bearing is the using it for ship 
main propulsion shaft line or thruster. In recent years fast 
development of engineering of polymer materials and their 
growing application to sea-water-lubricated slide bearings, has 
been observed. However the application of passive magnetic 
bearing of Halbach's configuration provides a comparably high 
stiffness of supporting, and simultaneously lack of wear-out 
phenomenon of interacting surfaces under operation. 

The magnetic bearing is a device which makes use of 
repelling forces between homopolar magnets simultaneously 
placed in front to each other and coaxially, which keep rotor 
against stator in the state of levitation. Multi-row arrangement 
of magnetic rings provides a higher stiffness of bearing unit and 
more stable work in the neigbourhood of point of operation [1]. 
The levitation means keeping the rotor supported without any 
mechanical contact between it and stator. The main advantage 
of such bearings is elimination of friction forces between 
interacting elements. In the commonly used rolling and sliding 
bearings friction phenomena as such is not eliminated but 
only reduced. Moreover, in the traditional bearing systems 


Fig. 1. Thruster based on a slip-ring motor 


additional systems responsible for lubricating, cooling and 
discharging wear products are necessary. As results from the 
presented characteristics the application of classical bearings to 
operation in the environments more aggressive to materials or 
ina high vacuum meets significant difficulties and limitations. 
Therefore it is recommended, wherever justified, to apply 
magnetic suspensions which do not require additional systems, 
can operate in chemically active environments, and eliminate 
friction, extending this way life time of traditional bearings. 
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Exchangable Blades 
in G.R.P. Material 


Fig. 2. Overall view of a Rim Drive device and the model intended for trials 


Magnetic bearings are split in two basic groups: active and 
passive ones. The first group covers active bearings in which 
magnetic force value is proportional to position of rotor centre 
against hole centre and this way inversely proportional to value 
of a gap between interacting parts. In active bearing system the 
following items can be always found [3]: 

e asensor intended for the measuring of rotor position within 
air gap, 

* acontroller which transforms information on rotor position 
into steering signals, 

* an amplifier which, on the basis of the steering signals, 
changes value of magnetic force by changing value of 
electric current flowing through active bearing coils in axial 
direction. 


As results from the presented design of magnetic bearing, 
the rotor position stabilization is realized by the system of feed- 
back between rotor position and magnetic force. 

Passive magnetic bearings belong to the second group. In 
design of such bearings there is no feed-back system between 
rotor position within air gap and value of magnetic force. In 
design of passive bearings permanent magnets are used to 
produce magnetic forces [4]. 

The main advantage of passive bearings is its simple 
design, very high efficiency (as such bearings do not absorb 
energy during operation and provide simultaneously magnetic 
levitation) and relatively low cost. Nevertheless the lack of 
a feed-back system does not ensure a stable position of rotor 
around a given working point and does not make it possible 
to control its position within air gap. Active bearings are to be 
used in any case where a decisive condition is to maintain high 
precision in rotor positioning. Moreover it is not possible to 
build a fully isolated magnetic suspension system on the basis 
of only passive magnetic bearings (acc. Earnshaw theorem) 
[5]. Therefore in any solution of magnetically suspended rotor 
at least one actuator is to be active, and the remaining can be 
passive ones [6]. 
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Hence for the bearing system of water propeller rotor it 
was proposed to apply two radial passive magnetic bearings 
and one axial active bearing. The bearing system is presented 
in Fig. 3. 


Electric motor stator 
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Magnets nets 
Electromagnet 
of axial 
bearing Electromagnet 
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bearing 
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bearing 


Fig. 3. Bearing system for water propeller rotor 


In this paper a mathematical model of radial passive 
magnetic bearing is presented. Evaluation of the model is 
necessary to elaborate design algorithms for passive magnetic 
bearings. 


PASSIVE MAGNETIC BEARING 


Two magnets or two sets of magnets are used to form 
a passive magnetic bearing [5]. In the simplest case two 
magnets repelling each other are used. One of them is rigidly 
connected with casing of rotary motor. It cannot displace. The 
other one is placed onto free rotor and can do, together with it, 
a complex motion being a sum of translational and rotational 
displacements. The moving magnet has six degrees of freedom 
and behaves as a magnetic compass needle which tries to be 
always lined up along external magnetic field lines (Fig. 4). 
As in the classical magnetic bearing there is a stator fixed in 
machine casing (motionless magnet) and a race placed in rotor 
(moving magnet). Between races of rolling bearing there are 
balls placed in bearing cage, and in slide bearing oil film appears 
between its races. In passive magnetic bearing a medium in 
which the bearing operates, is used. If a bearing operates in high 
vacuum conditions then between its races vacuum occurs, and 
a bearing immersed in water contains water between its races. 
This feature of magnetic bearing is one of its crucial merits. 
Obviously, apart from an environmental medium also magnetic 
field occurs between bearing stator and race. 


Moving magnet 


Motionless magnet 


Fig. 4. Arrangement of magnets in repelling system 


a) b) 


Ring sector 


of moving magnet 


It can be assumed that in passive magnetic bearing the 
permanent magnet which does not displace is a source of 
non-uniform magnetic field. And, the permanent magnet 
connected with the rotor lines up along lines of external 
magnetic field generated by the magnet fixed in rotary 
machine casing [7]. 

As the magnet fixed in rotor behaves like a magnetic 
needle it tries to rotate its north pole towards south pole of 
the motionless magnet and shift its south pole so as to place it 
close to north pole of the motionless magnet. Magnet placed 
in a non-uniform magnetic field does a complex motion and 
its trajectory depends on initial position of moving magnet 
against that motionless. If magnets point to each other with 
opposite poles then they will try to be mutually connected. After 
connection the two magnets will behave as one magnet which 
becomes a magnetic dipole. If an orientation of the magnets 
is different from that above mentioned, then the magnet will 
rotate until they point to each other with opposite poles and 
then connect mutually. The connection of the magnets means 
that they reach state of equilibrium, hence such mutual position 
of the two magnets appears stable [7]. 

Any rotary displacement of rotor within magnetic bearing 
and mutual connection of its magnets would lead to a failure 
of rotary machine and its bearing. Only the maintaining of 
magnets in the position in which repelling forces between 
magnets occur, makes bearing work correct. The mutual 
repelling process of the magnets is a state of instable 
equilibrium, i.e. the state in which the system maintains its 
stability but with no stability margin. 

The bearing intended for the transferring of radial loads 
consists of two magnets concentrically located (Fig. 5a). The 
bearing magnets are radially magnetized. The same magnetic 
pole takes place both on the outer circumference wall of the 
magnet fixed in machine rotor and on the inner circumference 
wall of the magnet fixed in machine casing. Such solution 
makes it possible to generate mutual repelling the magnets and 
this way magnetic levitation of the rotor to be suspended. As 
the manufacturing of uniform radially magnetized ring magnets 
is difficult, also bearings of rings built with segments, can be 
found. Such segment is a ring sector which is to be magnetized 
radially (Fig. 5b). Due to technological reasons it is not possible 
to produce an one-sided radial bearing. Radial bearings are 
always manufactured in a differential system. 


Ring sector 
of motionless 
magnet 


Fig. 5. Radial passive magnetic bearing 
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a) b) 


Outer circumferential wall 
Right side wall 


Fig. 6. Surface current constrained within side wall (6a), coordinate frame orientation in the magnet (6b) 


MODEL OF RADIAL PASSIVE BEARING 


Amper model was used for description of a model of magnetic interaction between two ring magnets. As assumed in the 
model in the magnet walls flow constrained currents being a source of magnetic field. Value of the constrained currents depends 
on the product of magnetization vector and that perpendicular to magnet wall [2]: 


Ken (1) 
where: 


m 
M - magnetization vector 


=> à 
n — vector perpendicular to magnet wall. 


Fig. 6 illustrates surface currents constrained within ring magnet. As results from the relation (1) the currents will occur in 
side walls only. They will not occur on the outer and inner circumferential walls of the magnet. 

In compliance with Lorentz principle [2] the magnetic interaction force depends on the vector product of the surface constrained 
current and the vector of magnetic induction which penetrates the moving magnet: 


F(x,y,z) = f (k x B) da (2) 


where: 
= 


Қ — surface constrained current flowing through a given sector of the moving magnet, da’ 
В — vector of magnetic induction which penetrates the moving magnet sector of the surface area da’ [in the point of the 
coordinates (x, y, z)]. 


Fig. 6a presents the ring magnet with depicted surface currents. In the ring magnets the currents appear both in the left and 
right wall (Fig. 6b). 
For radially magnetized ring magnet the magnetization vector is equal to [4]: 
M'-[M, Му 0]-[M'coso M'sinp 0] 


The normal vector can be expressed as: 


п= [0 0 1] 
The surface current in the left wall amounts to: 
_ ï j К 
К" = |М'соѕф' M'sinog' 0| = М'ѕіпфт — М'соѕф7 
0 0 1 B) 


К" = [M'sing —M'cosọ' 0] 
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The surface current in the right wall of moving magnet is described as follows: 
[1 j kà | _ 
К' =|M'cosg’ M’sing’ 0 | = —М'зїпф1 + M'cose’y 

0 0 —1 


(4) 
K', = [-M'sino' М'соѕф' 0] 


induction in the point of the coordinates (x, y, z), is equal to: 


— 


The vector product of the surface current flowing through the left wall of moving magnet and the magnetic 


ї j 
КХВ B=] M'sing’ 


> 


k 


0 = 
В,(х,у,2) By(x%y,z) B(x y,z) 


(5) 
В, (х,у, 2) М'соѕфї — В, (х, y,z)M'sing) + M (B, (х, y,z)sing + B, (x y, z)coso")k 


The surface area da’ of an elementary sector of magnet side wall can be estimated as follows (Fig.7): 


—M'cosg’ 


da’ = r'do'dr' (6) 
where: 
т? — the mean radius changeable within the interval г’ <т?< К? 
Фф” — the angle changeable within the interval 0 < ф < 2л 


I 
t 


== 
2 


Fig. 7. The sector considered оп the moving magnet side wall 


By taking into account the above given relations in the Lorentz force (2) the following value of the force on the left wall is 
obtained: 


Fi(xy,z) = – [| B, (x, y, 2M coso r'do'dr'1 — || B, (x, y, 2) M'sing'r'do'dr'] + 


+ || M (By (х, y, Z)sing + В, (х, у, z)coso")r'do'dr' К 


(7) 
In the same way the magnetic force value in the points located on the moving magnet right wall, can be found: 
Ер(х,у,2) = | | В, (х, у, 2) М'соѕфт'аф'агт | | В, (х, y, z) M'sing'r'do'dr'j + 


- [| w,6.v. 2)ѕіпф + В, (х, y, z)coso")r'do'dr'k 


(8) 
By using Biot-Savart principle magnetic induction in the point p of the coordinates (x,y,z) located in the space around 
motionless magnet (magnetic field source), can be determined. Points for which magnetic induction is calculated for further 
analysis are placed on the moving magnet. Biot-Savart principle for surface currents flowing through motionless magnets takes 
the following form [2]: 
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Brn) — Ho ( K(p)xP 
B(p) = f ——7- da 
" р 
" where: 
P - versor of vector directed opposite to the surface on which surface constrained current flows 
to the point for which value of magnetic induction has to be determined. 


E 
The vector P can be written as follows: 


P = |P|P 
The relation can be transformed into the form: 
5 Р 
Р = — 
|P] 


By accounting for Eq.(11) in Eq.(9) the relation for magnetic induction is obtained as follows: 


Bim — Ho f KONP 
B(p) = |р? da 


(9) 


(10) 


(11) 


(12) 


Value of the magnetic induction depends on the vector P, therefore its coordinates between the point located on the circular 
loop circumference and the spatial point p’ (x,y,z), should be estimated. The coordinates of the point p’ located on the circular 


loop circumference are as follows: 


Coordinates for the left wall, р? (x^, y, 77”) Coordinates for the right wall, р”, (б? y» т.) 


T 
Xj = rcos@ X р = rcosq 


у | = rsing у p = гѕіпф 
—h 


The coordinates of the vector Р for the left wall are equal to: 
Р, = [x — x y-y' z—z 
and, for the right wall, as follows: 
P, = [x 7 Xp yx. z — 7'y| 


where: 
x,y,z — coordinates of the point for which value of magnetic induction has to be determined. 


Finally the coordinates of the vector P take the form: 
= | h 
В = |х – гсоѕф y-rsino 2—5 
— | һ 
Рр = |х – гсоѕф y-rsinp z+ 2 


The module of the vector [Pi] takes the value: 


Р, 2 
В| = |х2 + y? + z? — 2г(хсозф + ysing) — zh + r? iL 


And, the vector IP. is equal to: 


P. 2 
|P | = [х2 + y? + z? — 2r(xcoso + ysing) + zh + r? + 


For the radially magnetized ring vector the magnetization vector amounts to: 


М = [М, M, 0] = [Мсоѕф Msing 0] 


The values of surface current іп the motionless magnet left and right walls are the following, respectively: 


Кү = [Msing —Mcoso 0] 


A 
| 


= [—Мѕіпф Мсоѕф 0] 
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(13) 


(14) 


The vector product of the surface current Кү and the vector Р, is equal to: 


K; x P; = |-Mcoso (z — =) —Msino (z - Э М(хсоѕф + уѕіпф — ?] 


In the similar way the vector product of the surface current E, and the vector P, can be determined: 


(15) 
ek h : h : 
Kp XP, —|Mcosp|zt7) Мѕіпф(2+ 5 —M(xcosqo + ysing — г) (16) 
For the analysis is used the magnet sector da limited by the arc of the angle d« and the vector increment dr. 
The surface area of the sector is equal to: 
da — rdodr 
аф 


‚= 


| 

1 
\ 
\ 


Fig. 8. The sector considered for the motionless magnet side wall 


The magnetic induction value in the point p of the coordinates (x, y, z) is equal to: 


E КүхР( 
В(х, у, 2) = H 


Kp xP, 
177 da + fe? da) 17 
MENS d 
By taking into (17) account the relations (13), (14), (15) and (16) the following magnetic induction value can be obtained: 
= Hu - E ré 
B(xy,z)- ax Cii + С1)1 + (C12 + C22)] + (Суз + C23)k} 


B(x, у, 7) = [Bx By B;] 


where: 
—Mcoso (z = =) r 
Си = ; dodr 
nè 
(x? + у? + z? — 2r(xcoso + ysing) — zh + г? + x) 
—Msinqo (z = >) r 
Сә = [| | p 
(x? + y? + 22 — 2r(xcoso + ysing) — zh + г? + al 

Суз = 


М(хсоѕф + уѕіпф — r)r 


3 dodr 


zy 
4 


(x? t y? + 22 — 2r(xcos@ + уѕіпф) — zh + r? + 
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Мсоѕф (z + >) r 
C21 = | 


3 dodr 


2\7 
(x? + у? + 22 — 2r(xcosq + уѕіпф) + zh + г? + Ty 


Msino (z + =) r 
C22 = | 


з dodr 


247 
(x? + у? + z? — 2r(xcos@ + уѕіпф) + zh + r? + = 


—M(xcos@ + уѕіпф — r)r 
=) 


= dodr 


(x? + у? + 22 — 2r(xcoso + уѕіпф) + zh + r? + 73] 


By solving the above given integrals and accounting for 
magnetic induction values in Lorentz force it is possible to 
determine lifting force of the bearing. 


SUMMARY 


- The elaborated mathematical model makes it possible to 
assess magnetic force value and direction. Such model is 
necessary for determination of maximum magnetic force 
generated by radial passive magnetic bearing. Moreover 
the model can be applied to assessing dynamic parameters 
of active magnetic suspension. 

- The presented approach leads to analytical solutions 
which make analyzing magnetic bearings by means of 
mathematical programs in a simple way, possible. The 
model has been adjusted to the Matlab-Simulink engineering 
software which enables analyzing dynamic processes. 

- The software based on the finite elements method (FEM), 
commonly used for solving static magnetic fields and 
magnetic forces do not make it possible to perform 
a comprehensive dynamic analysis of radial passive 
magnetic bearing. 

- Basing on the elaborated model these authors conduct, 
by means of Matlab software, design process of passive 
suspensions applicable to water propellers. The elaborated 
models and relevant design algorithms will be verified 
during experimental tests and analyses with the use of 
FEM-based programs (static analyses). 

- Elaboration of comprehensive mathematical models and 
then determination of their allowable simplification would 
make it possible to build algorithms capable of facilitating 
implementation of passive magnetic suspensions in 
practice. 

- Apart from the elaboration of passive models of radial 
magnetic bearing, it is intended to elaborate models 
of axial bearings. Such bearings will be applicable to 
a bearing system of water propeller rotor, based on another 
configuration of magnetic bearings. 
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Problems of welding in shipbuilding - an 
analytic-numerical assessment of the thermal 
cycle in HAZ with three dimensional heat source 
models in agreement with modelling rules 
Part Hl 


Non-linear analytic-numerical assessment 
of thermal cycle - examples 


Eugeniusz Ranatowski, Prof 


University of Technology and Life Science, Bydgoszcz 


ABSTRACT 


This part is continuation of PART II. Analytic solutions for the temperature distribution 

in HAZ — presented in the previous part of this article are transformed for computer 

calculation with used Mathcad programme. There are established algorithms in moving 

and stationary systems for thermal cycle calculating. Finally, a few analytical examples 
with use of C-I-N and D-E models are demonstrated. 


Keywords: welding; shipbuilding; welding in shipbuilding; thermal cycle; heat affected zone; heat source model 


ADAPTATION OF THE ANALYTICAL 
SOLUTIONS FOR NON-LINEAR 
COMPUTER CALCULATIONS 


The equations (37)+(42) of PART II of this article are 
algebraic form of linear heat flow solutions. 

In order to execute computer calculations with temperature 
dependent physical parameters: A, c. , p the above algebraic 
expressions must be transformed. For this purpose we will use 
calculations in Mathcad programme [1]. This programme is 
very useful for modelling and simulation of welding thermal 
process [2, 3]. 

Therefore the following assumptions were done: 

* heat source energy is being input to the metal during time 

At, not impulsively At—0. HS inputs are being summed up 

in points in distance Ax = v At. Considering this t’ = (j-1) 

At, (j = 1, 2, 3 ..n). 

* integrals were replaced by finished sums assuring sufficient 
exactness. Finally, the following computing expressions for 

linear heat flow solutions are obtained [4]: 


A. from Cylindrical-Involution-Normal heat source 
model 


* Stationary co-ordinates system : 


T(x, yo. zo t) 9 > if t « G- 1): At.0, 


j=l 


л.с,-(@—ехр(—К-5)) 4-a-k-(t-0G-DAt+1 


Ec (1a) 


"exp ; 
4-a:k(t-(J- DAO +1 


last 


УВ, "C-D; :exp[- a: (t - G- DAO]] 


jl 
* Moving co-ordinates system: 
T(x, y.z,th= Dif t « (-1)- At.0, 


ial 


q:k-K, 1 


n-c, (l—exp(-K-s)) 4-a-k-(t-(j—-Iat+1 


=ke ~(j-1)v- А) + y? Ib 
p k-(+vt-G-Dv-At) «y )| (0 


4-a-k(t- (j- А) 41 
last 


->B C, D, -expa -r (t- G- DADI} 
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B. from Double Ellipsoidal configuration of source 


* Stationary co-ordinates system 
n 


Tq: ¥os%Zast)= > iff « ( —1)-At,0, 


jl 


q.f,:3-43. At 
А (2:a-(t- (j - А) +а;). 
:(12-a-(t- (j -DAt) - b?) 


&,-v-G-DA — , 


4-a-(t- (j- DAD 2-2? 
3 
exp} — А + 
Yo 


e 
4-a- (t-G-DA)+ b 


(2a) 


(12-a-(t- (j -DAt) a2) 
(12-a-(t- (j -DAt - b?) 


(x, - v: G- А) 


4-a-(t- (j - DAD a 
-exp| — : 


+ 


уо 


с 
4a - G- DAO- b; 


+o Cj: Dj -exp a- (t7 G- DADI 


In order to execute computer calculations with temperature 
dependent physical parameters: i, €, P the above algebraic 
expressions must be transformed. 

Therefore the following assumptions were done: 

* heat source energy is being input during time At, not 
impulsively At — 0. Hs inputs are being summed up in 
points in distance Ax = v At. Considering this t = (j - 1) At 
(712,3,...S), 

* theintegrals are changed to an algorithm executing proper 
summing with physical parameters upon temperature 
change control, 

e asA(T), e, (T), p(T), a(T) values in defined increments are 
known - like shown in table 1, the matrices containing T and 
corresponding A(T), c (T), p(T), a(T) values are defined. 


With use of linear interpolation procedure, the continuous 
functions МТ), с (T), p(T), «(Т) were created and built-in inside 
calculation sheet. 

There are three main mini-procedures responsible for 
thermal cycle calculating. In the first of them initial values are 
presented (these values are specific for the given cycle)."Stab - 
time" parameter is estimated time needed for stabilisation of 
thermal field in moving co-ordinates system, “At” is duration 
time of every heat impulse being input. Therefore “$” gives the 
total number of heat impulses to be generated in order to obtain 
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* Moving co-ordinates system 
n 


T(x. y.z.t)= vif tt <(G-1)-At,0, 


Jel 


q-f,-3-43-At 


(12-a-(t—(j-lAt)+ a7) 
-(12-a-(t-(j-DAt)+ bF) 


E UN 
a 


(x*wt-G-DAD | 
4-a-(i— (j- DAt)4- a 
exp| – В 3 + 


у 


=== 
4-a-(t- G- DAD zb; 


q-£, 3-45. At' d 
ЕА 
-(12-a-(t-(j-IAt)+b?) 
(x * v(t - (j - DAOY 
4-a (t - G- DAD a! 


+ 


“exp — 


2 


y 


a QN 
4-a-(t- G- DA zb; 


last 


УВ, C; р, -exp az! (t - = DADR 


the summary thermal field as a result (this parameter is being 
used finally in the third mini-procedure). The estimation of 
thermal fields from several impulses is running with changeable 
values of a and X parameters according to Table 1. Parameter 
“Last” is used in the second mini-procedure which computes 
г, Т, ,,....r, Values — again with step by step X(T) values being 
modified. The final mini-procedure summarises thermal fields 
from several heat energy impulses using a proper formula 
specific for several HS model. 

The following algorithms (3) and (4) with initial computing 
parameters are usually used to perform calculations for various 
heat sources in moving co-ordinates system: 


A. from Cylindrical-Involution-Normal heat source 
model: 
Stab time — 19 sec At — 0.05 sec 
S = Stab_time/ At = 380 last = 0 

T<0 

Stab_ time 


for je1..S- 
At 

T, (x y. zt) 2|a << a(T) 

à < МТ) 


for ie1...last 


n « threshold 


[ cot(n-g) — 


X ЕЧ -0, d, 
root g 


ЕКЕН 
g 


П 


+ 


n < threshold 
[ cot(ng) – 


X ЕЧ -o 
root g 


П 


ЕКЕН 
5 


+@—-1).2 T mS 
E 8 " ер oo 8 = 
[ cot(n- g)— cot(n - g) — 
3 — 1: | ЫЙ a. 
if] root n ЕЯ Aes if| root A e 1) d мы 
E -n 
ЕКЕН ЕКЕН 
L g | L g 
> threshold > threshold 
threshold + (i - 1) - Ž otherwise threshold + (i 21) otherwise 
g g 
[T< Tif t«(j-D-At0,q- T«Tif t«(j-D-At,0, 
ket. 
чл q-f, -3-v3-At | 
— _ дА (I2-a-(t-(j-At)+a,")- 
"T =f dom] s TA —0.- 


[x * v-t- G-DAOY * y] 


* ex = |e eoe * 9 2 
"^ телок - G-DAO ais г 
lat | | 4-a-(t—(j-NAt)+—-a, 
-Š B; C, - D; expla -r -(t-(j-1)- At)] eal = 3 " 
i=l 2 
Т "— — — 
т ЕВА се и 
В. from D-E configuration of source: a) ке P 
T<0 q:£.-3-43-At 
for jel g — Stab_time + 3 
br jel.. E ji (12-a-(t- (j- DAD a, ^)- 
n-vn-—-c,- 2 
Т, (х,у,2,0) =|а <a(T) а :(12:a-(t- (j -DAt)- b, ^) 
= ACT) 
for ic1.las Gev-(t-G-DADY , 
4-a-(t- )А)+1.а2 
exp — 3 


2 


> __ 
das t7 Q- DAD b; 


+ 


Ув, C; Ej -expl -r «(E-0-D«A10] 


i=l 


For stationary system (x, = xtvt, y, = у, Z, = Z) the 


following algorithms (5) and () are presented: 
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C. from Cylindrical-Involution-Normal heat source 


model: 
Stab_time = 19 sec At = 0.05 sec 
S = Stab_time/ At = 380 last = 0 


T<0 
iujeio- Stab time 
At 
T, (x. y. zt) =|а < a(T) 
à «— МТ) 
for 1є1..Јаѕї 


n « threshold 


[cot(n:g) – 


2 А л ? 
X -|n+G-1)-—| -ao:a 
root [ pog d Mim 


ЕКЕН 
g 


П 


+ 


root 


root 


n « threshold 


[cot(n:g) – 


ZI -a,-0, 
g 


ЕКЕН 
g 


"n 


+@—-1).2 


|со(п-р)— 


X ЕЧ EREN 
8 


d 
| g 


> threshold 


«TY 


threshold+(i—1)-~ otherwise 
g 


А 4-0) 
i g " E А 
Є фы Я T«T if t<G-1)-At,0, 
. №. ЕЧЕ 
| 
"es : a q-f,-3-J3-At 
d - Я 
um! g dz A А (12:a-(t- (j- DAt)* a; )- 
TT: _—e . 
» threshold a — V -02-a-(t-G-DAO4 b) 
threshold + (1—1): otherwise : 
Р (,-v-G-DADY — | 
ТТ i . I. 
| 4-а-(@—-(0—1)Аї)+—-а, 
t<(j-1)-At,0,q: exp — , 3 РВ 
k— К,А! E EL 
n-—-(1—exp(-K, :5)) 4-a-(-G- DAD zb, 0 
a 
oe RIEN | КОСЕ, 
1+4-a-k(t—(j-l)-At , 
о d NN в) kg, [@-а-@-0-50ло+а,®); 
6 2 TXTU— с. 
-exp a Wty C G- DAD +y"] (2-a-(t- G7 DADO» b,”) 
1«-4a.-k(t- (j- DA 
last (x, - v: (G-DAOY 
D еб EPOR 
i=l ы а шуу" 
-exp[a r? «(t - (j- 1)- At)] -exp| — ; 
T — M 
D. from D-E configuration of source: 4-а-((— (j —NAt)+ 3 ` b. 
T<0 
last 
Я Stab_ time 2 3 
for je. == D CE; expli. -(t—(j-1)- At)] 
T, (x. y. zt) =|а < a(T) . . 
MT) The base procedure defined by algorithms (3) + () requires 
Я = cooperation with several sub-procedures such as Stab_time, 
or iél...last 
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last, calculation of roots г, , A(T), a(T). 


From essential point of view the sub-procedure Stab_time 
defines in time dimension the moment of time after which it 
estimated its value needed for stabilisation of thermal fields 
as follows: 

- Sub-procedure Stab time 
Stab time(x, y, Zz) = 
absolute < 20 
settled < 10 
stroke < 1 
a MOVE(x, Sf s absolute) 


«1.02 
T_MOVE(, у, z settled) 


while 
if 
settled «— settled — stroke 
x NEUE MOVE(x, у, z absolute) 
F- MOVE(x, y, Z, settled) 


otherwise 


>1.02 © 


settled < settled + stroke 


(settled + stroke) 
T MOVE(x, У, 7, absolute ) 
T_ MOVE(x, y, z settled + stroke) 
T. MOVE(x, у, z absolute) 


Stab_time is compound sub-procedure realising calculation 
of wanted value of stabilised time for optional point in moving 
coordinates system. Parameters “absolute” and “settled” are 
preliminary set up values of start and end of partition of time, in 
which follows search of time stabilisation with “stroke” step. 

The procedure is built in this way that last value of 
Stab_time does not depend on value “absolute” and “settled” 
but correct assumption of these values shorten time wanted for 
account of Stab_time. 

Furthermore in the first main mini-procedure also cooperation 
is required with sub-procedures for inside calculation of X(T), 
a(T). The discrete values of A(T), a(T) are known and shown 
in Table 1 the matrices containing T and corresponding A(T), 
a(T) values are defined experimentally. Than with use of linear 
interpolation procedure, continuous functions were created and 
built in inside calculation sheet as follows: 

- Sub-procedure A(T) 


(1) | (2) 
0 | 0.6285 
100 | 0.5866 
200 | 0.5447 


300 | 0.5028 


5% [0.5028] | (7) 
і < floor — 

400 | 0.4609 100 
500 | 0.419 i<14if T21500 
| 600 | 0.3771] |- (A 227 A 12) (8) 
700 |0.3487 — 

А icu 


^7 '800 [0.3268 Ал 
L (t-A Lu)* Aus 
T:-1..1500 


МТ): = 


900 | 0.3226 
1000 | 0.3268 
1100} 0.331 
1200 | 0.3352 
1300 | 0.3352 
1400 | 0.3352 
1500 | 0.3352 


- Sub-procedure о(Т) 


(1 (2) a(T): = 
Lee [ie floor Т ) 
100 | 0.16 100 
ba 
[2% [013] | т<тёв 
400 10.093] | |i€- 8«if T2 768 
500 |0.079 i<9<if T 2 800 
| 600 |0.062| | |i «- 10 if T2900 
C a “licll<if T2901 (9) 
800 0.042 i«12<if T21200 
| 900 | 0.055 (Ci.22-C int) 
901 | 0.062 €....-C,.,,) 
1200 | 0.062 n 
1300/0.062| | | (-С„)+С 
|1400]0.062] | T:=0...1500 
1500 | 0.062 


On Fig. 1 are presented discrete values of A(T), a(T) and 
determined by continuous functions with used sub-procedures 
(8) and (9) for a. A(T), b. o(T). 

We have high conformity of continuous functions and 
discrete value of A(T), a(T) from above-mentioned date on 
Fig. 1. 


А? 
оо 
(Т) 
03 5750 300 450 600 750 900 1050 1200 1350 1500 
a. дете 
© Discreet values of 2 (T) 
«+++ Course of function à (T) 
<2 
aa 
a (Т) 


0 150 300 450 600 750 900 1050 1200 1350 1500 
b. бү 

B Discreet values of ot (T) 
«+++ Course of function О. (Т) 


Fig. І. Values of a. A(T) and b. a(T) in agreement with tab. 1 
and continuous functions of for low carbon steel 0.1% C 
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The sub-procedure “LAST” realises calculation of value 
of parameter “last” defining correctness of calculation in 
analytical meaning but it decides indirectly about amount ofr, 
roots which is taken into consideration in each computational 
step. It is possible to assure required accuracy of calculation 
when value of “last” is sufficiently large. 

Sub-procedure “LAST” can find this value through 
analysis of moment of numeric convergence of neuralgic 
mathematical module in most adverse conditions at t = 0 and 
z = 0 as below: 

- Sub-procedure “LAST” 


i1 
_ || Totali үү + 
while} | Total; 
Last = + (Total, < 0.96: max(W)) | (10) 
iil 
iifi>15 
15 otherwise 
Pulse, = = cos(r, -Zo )+ 
(1- exp(- K -s))-— 
D (11) 
Ao - 
+ -sin(r, -Z,)-B,C,D, 
A-r; 
Total; = > Pulse, (12) 
i=l 


Realisation of sub-procedure “LAST” relies on search 
of sufficiently large value “i”. In order to set up condition 
of assumption converge of monotones growing series (from 
foundation for value 0.001) at necessity outreaching of 9% of 
maximum value “Total” in treated interval “i”. It catches on that 


minimal value of “i” and can’t be smaller than 15. On Fig. 2 it 
is presented course of value “Pulse” and “Total”. 


Fig. 2. Course of value “Pulse” and “Total” 
as function of growing value “i” 


Furthermore let’s notice that computational cycle is started 
with acceptance of value A(T = 0), «(Т = 0) but each next time 
step “At” change the values of A(T), a(T) in agreement with 
current value of temperature T. Generation of value “т is next 
computational step and for each time step, which is found in 
numeric way according to separate “sub-procedure г,” which is 
directly insert in algorithms (3)+() and simultaneously checks 
accuracy of calculation with used “threshold” parameter. The 
value of “threshold” parameter is established on very close 
zero but unzero. 

For example in Tab. 2 some values are presented of roots r, 
along with estimates of their accuracy for t= 0 and g = 1.2 cm, 
threshold = 0.0000001, a, = 0.02 W m?K'",, a, = 0.01. 

The final sub-procedure "Temperature — T" in algorithms 
(3) + () summarises thermal fields from several heat energy 
impulses using a proper formula specific for appropriated 
H-S model. 


EXAMPLES 


A few examples of welding cycles in stationary co-ordinates 
system are shown. C-I-N and D-E configuration of source are 
used for analysis. The results from C-I-N and D-E models 
are compared with Rosenthal-Rykalin solution along with 
experimental results? . Material parameters are accordance 
with Tab. 1 for low carbon steel. 


Tab. 1. A(T), c (T), p(T), and a(T) values in several temperatures for low carbon steel — 0.1%C 


рТ) (T) 
JK 'c™3 


50 


POLISH MARITIME RESEARCH, No 3/2010 


Tab. 2. The assessment of a few r,roots (11 of 60) for non-linear approach calculation at t = 0 


[D 


r, values 
0.198811490 


r, check 
5.153050313 - 10? 


l2|— 


2.633099949 
5.243573590 


3.597482134 - 10° 
6.028955113 · 10? 


AJU 


7.859042948 
10.475772574 


6.232056649 - 10° 
5.70471002 - 107 


ON | 


13.093007438 
15.710495155 


5.118749868 · 10? 
3.640161481 · 10° 


18.332817429 


2.0121055824 · 107 


О |оо| NI 


20.945850074 


1.066098321 - 10° 


23.563632979 


4.666276254 - 107? 


26.1814807 


3.501158403 - 10° 


Tab. 3. Results of estimate temperature 


Temperature [°C] 


No of 


1. experiment 2. C-I-N calcul. 


3. D-E calcul. 4. R-R calcul. 


example 


max max 10s 


20s max 10s max 10s 


Example 1. 

Main parameters: q = 2400W; g = 0.4 em; v = 0.75ст s” 
C-I-N: s = 0.2 cm; Kz = 15 ст; k = 12 cm? 

D-E: = 0.; a, = 0.5 cm; b, = 0.5 cm; c,= 3.0 ст; 

f = 1.4; а = 1.0 cm; b, = 0.5 cm; с = 3.0 cm — Fig.3. 


10 20 
t [sec] 

Fig. 3. Temperature change in points: x, = 0, y, = 0.1, z, = 0 
Symbols: 1 — experimental, 2 — C-I-N calculation, 3 — D-E calculation 
4 — Rykalin-Rosenthal solution (plate model with the sector line source) 


Example 2. 

Main parameters: q = 3300 №; g = 0.8 cm; v = 0.5 cm s 
C-I-N: s = 0. cm; Kz = 5 ст; k = 12 cm ? 

D-E: f.= 0.; а,= 0. cm; b, = 0. cm; с,= 4.2 cm; 

f = 1.4; a = 1.2 cm; b, = 0. cm; c, = 4.2 cm — Fig. 4. 


1500 | 


1000 


500 | 


1 


20 


10 15 
t [sec] 


Fig. 4. Temperature change in points: x, = 0, y, = 0.1, z, = 0 
Symbols: 1 — experimental, 2 — C-I-N calculation, 3 — D-E calculation 
4 — Rykalin-Rosenthal solution (plate model with the sector line source) 


Example 3. 

Main parameters: q = 12000W; g = 0.8cm; у = 1.1 cms" 
C-I-N: s = 0. cm; Kz = 2 cm; k = 8 cm? 

D-E: i 0.; a, = 0. ст; b,- 0.7 cm; c, = .2 em; 

f = 1.4; a = 1.4 cm; b, = 0.7 em; c, = .2 ст - Fig. 5. 
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Fig. 5. Temperature change in points: x, = 0, y, = 0.1, z, = 0 
Symbols: 1 — experimental, 2 — C-I-N calculation, 3 — D-E calculation 
4 — Rykalin-Rosenthal solution (plate model with the sector line source) 


For the comparison of account and experiment results, in 
Table 3, effects of estimated temperatures for examples 1 + 3 
are collected. 


Temperature is defined for: 
T,,,- maximum temperature, 


Т. — temperature after 10 5, 
Т, = temperature after 20 s. 


In the first example we can see distinctly correspondence 
of assessed temperature and experiment when we used 
analytic-numerical method and equations (1a), (2a) with use 
C-I-N and D-E heat source models respectively. For above 
example it take a stand of difference of order 100 - 300°C 
for maximum temperature. Results got run away from these 
issues is analytical assessment of run of temperature by used 
of Rosenthal-Rykalin (R-R) solutions. 

In the second example it takes a stand a similar situation 
but with certain difference in upper temperature. Highest 
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conformity with experiment in this range of temperature is 
for analytic-numerical solution and agreement with equation 
(1а) for C-I-N heat source model. Results of estimates 
temperature with used pure analytical R-R solutions definitely 
run away from analytic-numerical solution especially in upper 
temperatures where divergence in estimated T, amount to 
500°C. Results of third example are similar to second example 
in essential meaning. 

The special feature of characteristic of above examples 
1 +3 depends on heat source power and line energy of welding 
(1 — 3200 J cm“, 2 — 00 J cm" and 3 -10909 J cm’). 

Furthermore these issues indicate on utility of analytic- 
numerical solution with adopted C-I-N heat source model rather 
for simulation of welding process of high concentrated energy 
used but may also be used for the simulation of arc welding 
process similarly as D-E heat source model. It confirms also 
laser welding process simulation [7]. 


CONCLUSIONS 


In this work some extended consideration about analytic- 
numerical methods conforming has taken place. 
It is obvious that: 

* with an application of various heat source models one can 
obtain very effective temperature field solutions, 

• with appropriate algorithms, calculations are very attractive, 
effective and can be quickly executed on PC computers, 

* further impact should of course be put on still more detailed 
welding phenomenon analysis. The specificity of metal 
phase change and other complicated phenomena should 
be discovered and reflected in complex model in order to 
make more accurate and detailed analysis possible. 


The results of proposed methods were compared with 
experimental data and Rosenthal-Rykalin solution. The 
accuracy of D-E and C-I-N results having in mind experimental 
data was discussed and there's no doubt that the accuracy of old 
solutions (R-R) seems to be out of date for these examples. 

All this makes analytical solutions very competitive with 
numerical ones and makes them very useful in engineering 
practice. 
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NOMENCLATURE 

“absolute”, s, - preliminary assumption of upper value of 
time stabilised 

“settled”, s, - preliminary assumption of lower value of 
time stabilised 

“stroke” , s, - astep assumption on way of installation of 


time of stabilisation 
a temperature field in moving coordinate 
system, °C. 


T MOVE (х, y,z,t) - 
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Analysis of fracture toughness 
of structural timber 


Part I 
Theory and experimental tests 
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Polish Naval Academy in Gdynia 


ABSTRACT 


This paper presents fundamentals of fracture mechanics of anisotropic materials. Fracture toughness of 
anisotropic materials, e.g. structural timber, depends a. o. on state of stress, environment, temperature 
and changes due to ageing, therefore in such materials cracking process runs in various ways. Timber is 
characterized by eighteen coefficients which determine its fracture toughness in contrast to metals for which 
this number is as low as three. In this part of considerations a way of conducting the tests of specimens 
subjected to cracking as well results of the tests of natural and modified timber, are presented. 


Keywords: anisotropic material; anisotropic material cracking; CAD specimens; 
test stand; complex loading; stress distribution 


INTRODUCTION 


In the majority of structures operating under load the sudden 
cracking of elements occur in spite of that stresses in them do 
not exceed their allowable level. Such cracks can be caused 
both by material defects resulting from improper technological 
processes, and by external random overloading. 

This work is aimed at the determination of relations between 
fracture toughness of natural and modified wood and size 
of defects, in given loading conditions. Results of the tests 
should reveal an effect of modification of wood on its fracture 
toughness. 

In order to be able to determine structural material fracture 
toughness it is necessary to determine experimentally material 
constants, mainly the stress intensity factor and critical crack 
length. It was consisted in determining quantitative relations 
between resistance against the cracking in given loading 
conditions and size of defect. Fundamentals of fracture theory 
of orthotropic materials are presented. Ways are discussed of 
determining the critical load and possibilities in determining 
the structural timber fracture toughness in selected anatomic 
directions. 

The tests were conducted for the two crack propagation modes 
where stratification develops along fibres in crack containing 
specimen: tangential-longitudinal (TL), radial-longitudinal (RL). 
Specimens were subjected to load leading to crack (gap) opening 
(disruption). 


BASIC RELATIONS WHICH DESCRIBE 
CRACKING PROCESS OF ISOTROPIC 
MATERIALS 


Fracture mechanics is applied to assessing material 
resistance against loss of its coherence. As results from many 
literature sources concerning fracture of isotropic materials, 
crack-containing material can be loaded in three ways which 
may lead to crack opening (rupture), longitudinal shearing 
(slipping) or lateral shearing [5, 6]. The specified modes are 
called pure modes of loading the crack (gap) - containing 
element, marked I, II and III, respectively, (Fig. 1), [5, 6]. 


Fig. 1. Basic modes of cracking depending on a mode of loading 
applied to a body: I) Rupture (crack opening, pulling apart); 
II) Longitudinal shearing (slipping); HI) Lateral shearing 


Apart from the specified loading modes a few cases 


of complex loading which are combinations of the above 
mentioned cracking modes, can be met in reality. 
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For description of fracture toughness of specimens made 
of a given material the following parameters are used: critical 
value of the energy release coefficient G, stress intensity 
factor K, , and integral J. The parameter С is determined as 
the derivative of the potential energy released during cracking 
process, U, taken with respect to the crack length a [2]: 

G= oti (1) 
да 

The energy release coefficient corresponding to the loading 
modes I, II, Ш are marked G, Gi, С, respectively. As assumed, 
if material fracture toughness is only slightly dependent 
on crack length rate then in order the gap to be capable of 
developing in an uncontrollable way, the following condition 
is to be satisfied: 


G, E Gy (2) 


The condition (2) is assumed to be the crack propagation 
criterion, and the parameter G, which characterizes material 
fracture toughness, is determined experimentally[6]. The 
critical values of the energy release coefficients G,, Сб. Сү 
characterize material fracture toughness. The stress intensity 
factor K, is used as a measure of stress and displacement 
field in the vicinity of gap vertex. Its value depends on value 
of the external loading o,., crack length a as well as factor Y 
(a constant dependent on a crack opening mode and specimen 
geometry, [1, 3, 5]). 


K, =К (с, „а, Y) o=I, II, Ш (3) 


Та the subject-matter literature the stress intensity factor 
corresponding to the pure loading modes is marked К, К, 
K p respectively [2, 3]. To initiate cracking process, its value, 
like that of the stress intensity factor, should exceed the critical 
value K for a given material, in compliance with the following 


inequality [3,5]: 
K, 2 Ke (4) 


The relation (4), in contrast to the relation (2), concerns 
only the area around gap vertex. 

Between the energy release coefficient G, and stress 
intensity factor K, the following relation occurs [3, 5]: 


СЕЕ К (5) 
where: 


c, - coefficient depending on state of stress, loading mode 
and mechanical properties of a given material. 


The third parameter which determines resistance of 
materials to loss of coherence, is the so-called integral J.. 
The integral J , like the preceding parameters, is a measure of 
energy necessary to trigger crack developing. In contrast to 
the parameters G, and K , the parameter in question accounts 
for a significant plastic deformation occurring in structural 
materials [1, 3, 5]. 


BASIC RELATIONS WHICH DESCRIBE 
CRACKING PROCESS OF ANISOTROPIC 
MATERIALS 


Structural timber, an orthotropic material, is characterized 
by that it is subjected to cracking in variuos ways, depending 
On stress state, environment, temperature and ageing changes. 
The material is stratified and undergoes delamination under 
load. The delamination of such composite material results from 
many factors such as: assembling errors, mechanical failures, 
stress concentration in vicinity of defects nad cracks as well 
as changes in material structure during drying. 
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Timber is a non-homogenous, anisotropic (multi-directional) 
material, that makes decribing its strength properties difficult. 


Heart-wood 


White-wood Early-wood 


Late-wood 


Transverse 
cross-section 


Longitudinal - radial 
cross-section 


Longitudinal - tangential 
cross-section 


Annual ring Longitudinal direction 


Radial direction 
R, xi 


Fig. 2. A stump of wood with marked anatomic directions 


As arrangement of rings in a given volume is differentiated, 
three principal directions (three planes of symmetry) are 
specified in wood: radial — R (x,), tangential — T (x,) and 
longitudinal one - L (x,), Fig. 2. If a wood specimen is cut 
sufficiently far from the centre of wood stump, so as to be able 
to consider curvature of rings negligible, then its properties 
can be deemed orthotropic. Wood considered as an orthotropic 
material have three characteristic directions strictly associated 
with its anatomic structure [4]. 

For orthotropic materials six basic modes of crack 
propagation are distinguished (Fig. 3). The first letter of 
the symbols stands for the direction normal to plane of 
crack, the second - direction of its propagation. Hence the 
tangential - longitudinal mode of propagation TL, tangential- 
longitudinal TR, radial-tangential RT, radial-longitudinal RL, 
longitudinal-tangential LT and longitudinal-radial LR, are 
distinguished [4]. 

For three modes of loading eighteen possible modes of crack 
propagation are assumed. Two modes are distinguished where 
delamination develops along fibres: tangential — longitudinal 
— TL and radial-longitudinal — RL, as well as four ones in lateral 
direction. In practice, structures are subjected to complex loads, 
e.g. ТЇП or WHI. 

For assessment of fracture toughness of orthotropic 
materials linear fracture mechanics is used as any standard 
method is still lacking. In fracture toughness testing, critical 
values of the coefficients G,» С, Сү and factors K,, Kio 
K n have to be determined for every propagation mode. Wood 
is characterized by eighteen coefficients which determine its 


Fig. 3. Crack propagation modes in orthotropic material 


fracture toughness in contrast to metals for which the number 

is reduced to three. 

An important problem of fracture toughness of materials 
is determination of the critical load Р. The critical load is 
defined as such value of force under which loss of specimen 
material coherence occurs. Its determination is complex. For 
determining P. value a few methods are used out of which the 
following can be distinguished [6, 7, 8]: 

- the method NL in which it is assumed that initiation of 
cracking process is commenced in the instant when a non- 
linearity of the P-6 (force- gap opening) diagram appears 
(Fig. 4); 


a 


Pas) 


Peal) 


Pega) 


Load P 
Signal EA 


Displacement б 


Fig. 4. Methods for determining the critical load P, 


- the method Р(ӧ) in which the point of intersection of the 
diagram P(6) and the straight line whose slope angle tangent 
is smaller by 5 % than that of the straight line which contains 
the linear part of the diagram P(6), is assumed the beginning 
of cracking process (Fig. 4); 

- the method based on non-linear and linear approximation 
of the diagram Р(8), (Fig. 5) 

- the method based on acoustic emission (EA) [6]. 


TESTS ON STRUCTURAL TIMBER 
CRACKING 
Geometry of specimens 


Tests of fracture toughness of composite materials can be 
performed with the use of specimens of different geometrical 


Pear) 


Load P 


Displacement б 


Fig. 5. Method for determining the critical load P , 
based on approximation of the diagram Р(д) 


forms (Fig. 6) [6]. In structural timber testing the double 
cantilever beam specimens (DCB) were used. 

The critical energy release coefficient G,, [2] for a given 
specimen is described by the relation: 


2 
6.22 (6) 
2 B да 
where: 
P, — critical load 
C — specimen flexibility 
B — thickness (breadth) 
a — crack (gap) length. 


The critical energy release coefficient С | for an assumed type 
of specimen is determined from the following formula [2]: 


_ э Ред 
E 0 Bea 


G (7) 
where: 
ӧ — crack (gap) opening. 

The corrected crack length equal to a + |A| was determined 
on the basis ofthe beam theory in which experiment conditions 
and influence of real geometry have been accounted for. The 
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a) с) 
Р Р 


m P 


d) e) P 


P 


Fig. 6. Specimens for determining fracture toughness of wood elements loaded in compliance with the mode I: a) double cantilever beam specimen, 
b) one-sided -notched specimen under three- point bending, c) one-sided-notched specimen under tension, 
d) double-sided- notched specimen under tension, e) specimen for splitting test 


unknown quantity A was experimentally determined from 
the diagram C!%(a). The relation was built graphically by 
measuring the flexibility C for a few different values of the 
crack length a. The straight line approximating the diagram 
С'З(а) determines the quantity A on the horizontal axis level 
(Fig. 7). The way of determining the quantity A is descibed in 
detail in the work [7]. 


УА 


om a 
A 


Fig. 7. Method for determining the linear correction A 
to the beam theory[3] 


For the double cantilever beam specimens the critical 
energy release coefficient G, is determined from the formula 
as follows [7]: 

3 P. 


* ^72 B.-(a«|Ap (8) 


The specimen flexibility C is determined by using Berry 
method, 1.e. by means of the relation: 


C=K-a (9) 


where: 

n —acoefficient experimentally determined on the basis of 
the diagram In(C) versus In(a) 

K — numerical value determined from the point of intersection 
of the diagram with the ordinate axis. 


The diagram In(C) versus In(a) can be built by measuring the 
flexibility C for different values of the crack length a (Fig. 8). 
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The coefficient n is equivalent to its slope angle tangent. The 
critical energy release coefficient С, is determined by using 
the following relation [8]: 


‚Р. 
(p е (10) 
2 Ba 
lg C 
Slope n 


lga 
Fig. 8. Method for determining the coefficient n 


Characteristics of the tested material and 
description of the test running 


A natural and modified pine wood was used in the tests. 
Modification of wood by means of methyl polymethacrilate 
improves its properties. The detailed decription of run of the 
modification process by using methyl polymethacrilate as well 
as properties ofthe modified wood are given in the work [4]. In 
Tab. 1 апа 2 are exemplified values of material constants as well 
as strength properties of the natural and modified wood [4]. 

The fracture toughness tests of natural and modified wood 
were conducted with the use of the universal electrohydraulic 
testing machine MTS-81012. The tested specimens were made 
of the natural wood of abt. 1296 humidity and modified one 
of abt. 8% humidity. The tests were performed in the room 
temperature of abt. 20? C. The DCB specimens were used in 
the tests (Fig. 9). They were made of the same batch of wood 
on the basis of which their strength properties were earlier 
determined [4]. For the specimens the square cross-section of 
20 mm x 20 mm sides, was assumed. The gap was cut by using 


a saw of 1.8 mm in thickness and 70 mm in length, and then it 
was deepened up to 80 mm by using a thin blade. 


Tab. 1. Values of the material constants of the natural and modified wood 


Material constant 
E, = E, [GPa] 
E, = E, [GPa] 
E, = E, [GPa] 

G,, = G,, [GPa] 
С=О, [GPa] 
Gpr = G,, [GPa] 


Notation: the symbol K0.0+K0.56 stands for amount of methyl 
methacrilate [kg] per 1 kg of dry wood 


Tab. 2. Strength limits for natural and modified wood under tension, 
compression and shear load, respectively 


Kind of К, R. Е Au R, 

material | [MPa] | [MPa] | [MPa] | [MPa] | [MPa] 
K0.0 95 55 4.5 85 | 22.10 
K0.56 118 98 9 32 30.20 


Notation: R,, - tensile strength along fibres; R, - compression 
strength along fibres; R,,, - tensile strength across fibres; К, - 
compression strength across fibres; R, - shear strength 


Fig. 9. DCB specimen used in the tests, of the dimensions: 
1= 200 mm, a = 80 mm, c = 10 mm, B = 2h = 20 mm 


For description of wood fracture toughness the energy 
release coefficients G were used in view of anisotrophy of 
the tested material. The tests were conducted for two modes of 
crack propagation in wood, i.e. for the tangential — longitudinal 


direction TL and the radial — longitudinal direction RL. They 
were performed for the gap length a = 80 mm. The specimens 
complied with the provisions of the PN-EN 408 standard. The 
pine wood was suitably dried and prepared to testing. For the 
specimens 200 mm length was selected, i.e. ten times greater 
than the greater dimension of their cross-section which was 
on average equal to 20 mm in the case of the modes RL and 
TL. Load was applied to the specimens at the constant rate 
of 0.03 mm/s. Special holding devices (chucks) had to be 
manufactured to make performing the tests with the use of 
the MTS machine possible. The tested material was glued to 
the chucks by means of the ,,PATTEX'", a universal glue for 
wood. 

Each of the specimens was measured with the accuracy of 
0,01 mm and marked in line with its direction of orientation. On 
each specimen measuring length was made with the accuracy up 
to 0,01 mm, in order to analyze fracture toughness coefficient. 
During the tests values of the load P and displacement 8 were 
recorded. Value of the critical load was determined by means 
of the method P(ò) (Fig. 4). The load P and displacement 6 
was recorded by using the measuring system consisted of 
a computer and A/D card. 

For determining the critical fracture toughness coefficient 
G» it is necessary to know the coefficient n whose value for 
a given kind of wood is determined experimentally on the 
basis of the approximation of function of flexibility versus gap 
length. As in the available literature sources the n-coefficient 
value for Douglas fir equal to n = 2.80 is found and pine wood 
properties are close to those of Douglas fir, hence for further 
calculations n = 2.80 was assumed for the pine wood [9]. On 
the basis of the relation (11) were determined critical values of 
the energy release coefficient С, of the natural and modified 
wood specimens for the directions RL and TL. 

Based on the assumed value of the coefficient n = 2.80, the 
critical coefficient G,, of the tested and modified wood for the 
modes TL and RL can be determined from the relation: 


(1 2,80-P,-8 
a g B-a 


The tests of natural and modified wood fracture toughness 
for the mode TL demonstrated that the crack initiation was 
triggered off before the load applied to the specimens obtained 
its maximum value. The determined critical load was equal 
on average to 80% of the maximum load. The initial cracking 
run was mild, then passing to violent. As observed, in the 
natural wood specimens occurred a single crack which was 
next developing until the maximum stresses were exceeded 
(Fig. 10). In the case of the mode TL the crack propagated 
almost parallel to direction of the fibres. 

In the case of the mode RL the crack was initiated, like 
in the preceding case, already before the element reached 
its maximum load. The determined average load was equal 
to 85% of the maximum one. In the mode RL the crack was 
developing in a stable way both in the specimens of natural 
and modified wood, as shown in Fig. 10. The specimens 
under tension in the mode RL, cracked along the fibres. The 
crack propagation direction was in line with the anatomic 
direction of the fibres. The crack propagated along the 
hardwood-softwood boundary. In the mode TL, like in the 
mode RL, the crack propagated up to the moment of splitting 
the specimens, however in the specimens of modified wood 
it was propagating slowly. 

Results of the tests of the specimens made of natural and 
modified wood, loaded in the way I of the modes TL and RL, 
will be presented in the next part of the paper. 
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RESULTS OF FRACTURE TOUGHNESS initiates cracking process. The method is capable of 


TESTS OF STRUCTURAL TIMBER determining, in a simple and sufficiently accurate way, the 
instant of initiating the cracking process. As results from 
The cracks in wood specimens have developed in different the performed tests, the ratio of the load in question and 
ways depending on a propagating mode. For the mode TL, the maximum load was equal to about 0,80 for the mode 
value of the applied axial force was increasing up to the critical TL and about 0.85 for the mode RL. 
load P. In this mode, single cracks appeared in the specimens; - The fracture toughness of natural and modified pine wood 
they were developing until the maximum stress values were was determined experimentally by means of critical values 
reached. As shown in Fig. 10, the cracks propagated parallel to of the energy release coefficients G, complying with the 
direction of fibres. It was observed that cracks were developing relation (11), to be presented in the next part of the work. 
in the same way іп the tested specimens for the same loads - In the tests were used the beam specimens which make it 
and constant time intervals. The 20 mm long fracture was possible to model real phenomena of cracking the wood 
formed. The identical crack development process in the elements containing defects. The critical coefficients 
tested specimens results from the fact that they were cut from G,, were determined with the use of the beam theory 
neighbouring areas of the tested material. accounting for material anisotrophy. The critical values of 


the coefficients G, , determined for the splitting of natural 
and modified pine wood, were practically identical for both 
analyzed crack propagation modes: RL and TL. 


NOMENCLATURE 


crack length 

specimen thickness 

- energy release coefficient 
- stress intensity factor 

- longitudinal direction 
radial direction 

- tangent direction 

- specimen length 

- load 

potential energy released during cracking 
crack opening. 
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The Ship Handling Research and Training Centre at Ilawa is owned by the Foundation for Safety of Navigation 
and Environment Protection, which is a joint venture between the Gdynia Maritime University, the Gdansk University 
of Technology and the City of Ilawa. 


Two main fields of activity of the Foundation are: 


> Training on ship handling. Since 1980 more than 2500 ship masters and pilots from 35 countries were trained at 
Iława Centre. The Foundation for Safety of Navigation and Environment Protection, being non-profit organisation 
is reinvesting all spare funds in new facilities and each year to the existing facilities new models and new training 
areas were added. Existing training models each year are also modernised, that's why at present the Centre represents 
a modern facility perfectly capable to perform training on ship handling of shipmasters, pilots and tug masters. 


Research on ship's manoeuvrability. Many experimental and theoretical research programmes covering different 
problems of manoeuvrability (including human effect, harbour and waterway design) are successfully realised at 
the Centre. 

The Foundation possesses ISO 9001 quality certificate. 


Why training on ship handling? 


The safe handling of ships depends on many factors - on ship's manoeuvring characteristics, human factor (operator 
experience and skill, his behaviour in stressed situation, etc.), actual environmental conditions, and degree of water 
area restriction. 


Results of analysis of CRG (collisions, rammings and groundings) casualties show that in one third of all the 
human error is involved, and the same amount of CRG casualties is attributed to the poor controllability of ships. 
Training on ship handling is largely recommended by IMO as one of the most effective method for improving the 
safety at sea. The goal of the above training is to gain theoretical and practical knowledge on ship handling in a wide 
number of different situations met in practice at sea. 


For further information please contact: 
The Foundation for Safety of Navigation and Environment Protection 


Head office: Ship Handling Centre: 
36, Chrzanowskiego Street 14-200 ILAWA-KAMIONKA, POLAND 
80-278 GDANSK, POLAND tel./fax: +48 (0) 89 648 74 90 
tel./fax: +48 (0) 58 341 59 19 e-mail: office@ilawashiphandling.com.pl 

e-mail: office@portilawa.com 
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GDANSK UNIVERSITY OF TECHNOLOGY 


is the oldest and largest scientific and technological academic institution in 
the Pomeranian region. The history of Gdansk University of Technology is 
marked by two basic dates, namely: October 6, 1904 and May 24, 1945. 


The first date is connected with the beginning of the technical education at 
academic level in Gdansk. The second date is connected with establishing of 
Gdansk University of Technology, Polish state academic university. Gdansk 
University of Technology employ 2,500 staff, 1,200 whom are academics. 
The number of students approximates 20,000, most of them studying 


full-time. Their career choices vary from Architecture to Business and С 0) а f) 5 K U f) ive П jy 
9f Teennology 


Management, from Mathematics and Computer Science to Biotechnology 
and Environmental Engineering, from Applied Chemistry to Geodesics and 
Transport, from Ocean Engineering to Mechanical Engineering and Ship 
Technology, from Civil Engineering to Telecommunication, Electrical and 
Control Engineering. Their life goals, however, are much the same - to 
meet the challenge of the changing world. The educational opportunities адаг 
offered by our faculties are much wider than those of other Polish Technical + 
universities, and the scientific research areas include all of 21st Century Е It f O woio 

technology. We are one of the best schools in Poland and one of the best d G u O G 5 а f. + 
known schools in Europe — one that educates specialists excelling in the 


programming technology and computer methods used in solving complicated r : 2 ( 
scientific, engineering, organizational and economic problems. e [| € | | [| t5 IS Id f g w 
THE FACULTY OF OCEAN ENGINEERING 
AND SHIP TECHNOLOGY 


The Faculty of Ocean Engineering and Ship Technology (FOEST) as 
the only faculty in Poland since the beginning of 1945 has continuously 
been educating engineers and doctors in the field of Naval Architecture 

and Marine Technology. 


The educational and training activities of FOEST are supported by 
cooperation with Polish and foreign universities, membership in different 
international organizations and associations, as well as participation in 
scientific conferences and symposia. Hosting young scientists and students 

from different countries is also a usual practice in FOEST. 


The activities of Faculty departments are related to: mechanics and 
strength of structures, hydromechanics, manufacturing, materials and system 
quality, power plants, equipment and systems of automatic control, mostly 

in shipbuilding, marine engineering and energetic systems. 


FOEST is a member of such organizations like WEGEMT; The 
Association of Polish Maritime Industries and the co-operation between 
Nordic Maritime Universities and Det Norske Veritas. The intensive teaching 
is complemented and supported by extensive research activities, the core 
of which is performed in close collaboration between FOEST staff and 
industry. We take great care to ensure that the applied research meet both 
the long term and short term needs of Polish maritime industry. FOEST 
collaborates with almost all Polish shipyards. Close links are maintained 
with other research organizations and research institutions supporting the 
Polish maritime industry, such as Ship Design and Research Centre and 
Polish Register of Shipping, where several members of the Faculty are also 

members of the Technical Board. 


The Faculty of Ocean Engineering and Ship Technology is a unique 
academic structure, which possesses numerous highly qualified and 
experienced staff in all above mentioned specific research areas. Moreover, 
the staff is used to effective co-operation and exchange of ideas between 
specialists of different detailed areas. This enables a more integrated and 
comprehensive treatment of research and practical problems encountered 
in such a complicated field of activity as naval architecture, shipbuilding 

and marine engineering. 


The staff of the Faculty has strong international links worldwide, being 
members or cooperating with international organizations like International 
Maritime Organization IMO, International Towing Tank Conference ITTC, 
International Ship and Offshore Structures Congress ISSC, International 
Conference on Practical Design of Ship and other floating Structures PRADS 

just to name a few. 


GDANSK UNIVERSITY OF TECHNOLOGY 
Faculty of Ocean Engineering and Ship Technology 
11/12 Narutowicza Street, 80-952 Gdansk, Poland 
Tel (+48) 58 347 1548 ; Fax (+48) 58 341 4712 
e-mail: sekoce@pg.gda.pl 
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